TABLE OF CONTENTS
- ABINIT/m_ebands
- m_ebands/ebands_apply_scissors
- m_ebands/ebands_calc_nelect
- m_ebands/ebands_chop
- m_ebands/ebands_copy
- m_ebands/ebands_downsample
- m_ebands/ebands_enclose_degbands
- m_ebands/ebands_expandk
- m_ebands/ebands_free
- m_ebands/ebands_from_dtset
- m_ebands/ebands_from_hdr
- m_ebands/ebands_get_bandenergy
- m_ebands/ebands_get_bands_e0
- m_ebands/ebands_get_bands_from_erange
- m_ebands/ebands_get_carriers
- m_ebands/ebands_get_edos
- m_ebands/ebands_get_edos_matrix_elements
- m_ebands/ebands_get_erange
- m_ebands/ebands_get_gaps
- m_ebands/ebands_get_jdos
- m_ebands/ebands_get_minmax
- m_ebands/ebands_get_muT_with_fd
- m_ebands/ebands_get_occupied
- m_ebands/ebands_get_valence_idx
- m_ebands/ebands_has_metal_scheme
- m_ebands/ebands_init
- m_ebands/ebands_interp_kmesh
- m_ebands/ebands_interp_kpath
- m_ebands/ebands_interpolate_kpath
- m_ebands/ebands_move_alloc
- m_ebands/ebands_ncwrite
- m_ebands/ebands_ncwrite_path
- m_ebands/ebands_nelect_per_spin
- m_ebands/ebands_print
- m_ebands/ebands_print_gaps
- m_ebands/ebands_prtbltztrp
- m_ebands/ebands_prtbltztrp_tau_out
- m_ebands/ebands_report_gap
- m_ebands/ebands_set_extrael
- m_ebands/ebands_set_fermie
- m_ebands/ebands_set_scheme
- m_ebands/ebands_sort
- m_ebands/ebands_update_occ
- m_ebands/ebands_vcbm_range_from_gaps
- m_ebands/ebands_write
- m_ebands/ebands_write_bxsf
- m_ebands/ebands_write_gnuplot
- m_ebands/ebands_write_nesting
- m_ebands/ebands_write_xmgrace
- m_ebands/edos_free
- m_ebands/edos_get_carriers
- m_ebands/edos_ncwrite
- m_ebands/edos_print
- m_ebands/edos_t
- m_ebands/edos_write
- m_ebands/gaps_free
- m_ebands/gaps_print
- m_ebands/gaps_t
- m_ebands/get_eneocc_vect
- m_ebands/get_gaps_
- m_ebands/jdos_free
- m_ebands/jdos_ncwrite
- m_ebands/jdos_t
- m_ebands/klinterp_eval_bsd
- m_ebands/klinterp_free
- m_ebands/klinterp_new
- m_ebands/klinterp_t
- m_ebands/pack_eneocc
- m_ebands/put_eneocc_vect
- m_ebands/unpack_eneocc
ABINIT/m_ebands [ Modules ]
NAME
m_ebands
FUNCTION
This module contains utilities to analyze and retrieve information from the ebands_t.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, MJV, BXu) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
1) Remove npwarr, istwfk. 2) Use 3d arrays for ebands%nband 3) Solve issue with Hdr dependency
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 MODULE m_ebands 28 29 use defs_basis 30 use m_errors 31 use m_abicore 32 use m_xmpi 33 use m_htetra 34 use m_nctk 35 use netcdf 36 use m_hdr 37 use m_krank 38 use m_skw 39 use m_kpts 40 use m_sort 41 use m_dtset 42 use m_yaml 43 44 use defs_datatypes, only : ebands_t 45 use m_copy, only : alloc_copy 46 use m_io_tools, only : file_exists, open_file 47 use m_time, only : cwtime, cwtime_report 48 use m_fstrings, only : tolower, itoa, sjoin, ftoa, ltoa, ktoa, strcat, basename, replace 49 use m_numeric_tools, only : arth, imin_loc, imax_loc, bisect, stats_t, stats_eval, simpson, simpson_int, wrap2_zero_one, & 50 isdiagmat, get_diag, interpol3d_0d, interpol3d_indices 51 use m_special_funcs, only : gaussian 52 use m_geometry, only : normv 53 use m_cgtools, only : set_istwfk 54 use m_pptools, only : printbxsf 55 use m_occ, only : getnel, newocc, occ_fd 56 use m_nesting, only : mknesting 57 use m_crystal, only : crystal_t 58 use m_bz_mesh, only : isamek, kpath_t, kpath_new 59 use m_fftcore, only : get_kg 60 61 implicit none 62 63 private 64 65 ! Helper functions 66 public :: pack_eneocc ! Helper function for reshaping (energies|occupancies|derivate of occupancies). 67 public :: get_eneocc_vect ! Reshape (ene|occ|docdde) returning a matrix instead of a vector. 68 public :: put_eneocc_vect ! Put (ene|occ|doccde) in vectorial form into the data type doing a reshape. 69 public :: unpack_eneocc ! Helper function for reshaping (energies|occupancies|derivate of occupancies). 70 71 ! Ebands methods 72 public :: ebands_init ! Main creation method. 73 public :: ebands_from_hdr ! Init object from the abinit header. 74 public :: ebands_from_dtset ! Init object from the abinit dataset. 75 public :: ebands_free ! Destruction method. 76 public :: ebands_copy ! Copy of the ebands_t. 77 public :: ebands_move_alloc ! Transfer allocation. 78 public :: ebands_print ! Printout basic info on the data type. 79 public :: ebands_get_bandenergy ! Returns the band energy of the system. 80 public :: ebands_get_valence_idx ! Gives the index of the (valence|bands at E_f). 81 public :: ebands_get_bands_from_erange ! Return the indices of the mix and max band within an energy window. 82 public :: ebands_vcbm_range_from_gaps ! Find band and energy range for states close to the CBM/VBM given input energies. 83 public :: ebands_apply_scissors ! Apply scissors operator (no k-dependency) 84 public :: ebands_get_occupied ! Returns band indices after wich occupations are less than an input value. 85 public :: ebands_enclose_degbands ! Adjust band indices such that all degenerate states are treated. 86 public :: ebands_get_bands_e0 ! Find min/max band indices crossing energy e0 87 public :: ebands_get_erange ! Compute the minimum and maximum energy enclosing a list of states. 88 public :: ebands_nelect_per_spin ! Returns number of electrons per spin channel 89 public :: ebands_get_minmax ! Returns min and Max value of (eig|occ|doccde). 90 public :: ebands_has_metal_scheme ! .True. if metallic occupation scheme is used. 91 public :: ebands_write_bxsf ! Write 3D energies for Fermi surface visualization (XSF format) 92 public :: ebands_update_occ ! Update the occupation numbers. 93 public :: ebands_set_scheme ! Set the occupation scheme. 94 public :: ebands_set_fermie ! Change the fermi level (assume metallic scheme). 95 public :: ebands_set_extrael ! Add extrael to initial number of electrons to simulate e/h doping. 96 ! (assume metallic scheme). 97 public :: ebands_get_muT_with_fd ! Change the number of electrons (assume metallic scheme). 98 public :: ebands_calc_nelect ! Compute nelect from Fermi level and Temperature. 99 public :: ebands_report_gap ! Print info on the fundamental and direct gap. 100 public :: ebands_ncwrite ! Write object to NETCDF file (use ncid) 101 public :: ebands_ncwrite_path ! Dump the object into NETCDF file (use filepath) 102 public :: ebands_write_nesting ! Calculate the nesting function and output data to file. 103 public :: ebands_expandk ! Build a new ebands_t in the full BZ. 104 public :: ebands_downsample ! Build a new ebands_t with a downsampled IBZ. 105 public :: ebands_chop ! Build a new ebands_t with selected nbands. 106 public :: ebands_get_edos ! Compute e-DOS from band structure. 107 public :: ebands_get_jdos ! Compute electron joint-DOS from band structure. 108 109 public :: ebands_get_edos_matrix_elements ! Compute e-DOS and other DOS-like quantities involving 110 ! vectorial or tensorial matrix elements. 111 112 public :: ebands_interp_kmesh ! Use SWK to interpolate energies on a k-mesh. 113 public :: ebands_interp_kpath ! Interpolate energies on a k-path. 114 public :: ebands_interpolate_kpath 115 116 public :: ebands_prtbltztrp ! Output files for BoltzTraP code. 117 public :: ebands_prtbltztrp_tau_out ! Output files for BoltzTraP code, 118 public :: ebands_write ! Driver routine to write bands in different txt formats. 119 public :: ebands_get_carriers ! Compute carrier concentration from input Fermi level and list of Temperatures.
m_ebands/ebands_apply_scissors [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_apply_scissors
FUNCTION
Apply a scissor operator of amplitude scissor_energy.
INPUTS
scissor_energy=The energy shift
OUTPUT
SIDE EFFECT ebands<ebands_t>=The following quantities are modified: %eig(mband,nkpt,nsppol)=The band structure after the application of the scissor operator %fermi_energy
SOURCE
1722 subroutine ebands_apply_scissors(ebands, scissor_energy) 1723 1724 !Arguments ------------------------------------ 1725 !scalars 1726 real(dp),intent(in) :: scissor_energy 1727 class(ebands_t),intent(inout) :: ebands 1728 1729 !Local variables------------------------------- 1730 integer :: ikpt,spin,ival,nband_k 1731 real(dp) :: spinmagntarget_ 1732 character(len=500) :: msg 1733 !arrays 1734 integer :: val_idx(ebands%nkpt,ebands%nsppol) 1735 ! ************************************************************************* 1736 1737 ! Get the valence band index for each k and spin 1738 val_idx(:,:) = ebands_get_valence_idx(ebands) 1739 1740 do spin=1,ebands%nsppol 1741 if (any(val_idx(:, spin) /= val_idx(1, spin))) then 1742 write(msg,'(a,i0,a)')& 1743 'Trying to apply a scissor operator on a metallic band structure for spin: ',spin,& 1744 'Assuming you know what you are doing, continuing anyway! ' 1745 ABI_COMMENT(msg) 1746 !Likely newocc will stop, unless the system is semimetallic ? 1747 end if 1748 end do 1749 1750 ! Apply the scissor 1751 do spin=1,ebands%nsppol 1752 do ikpt=1,ebands%nkpt 1753 nband_k = ebands%nband(ikpt+(spin-1)*ebands%nkpt) 1754 ival = val_idx(ikpt,spin) 1755 1756 if (nband_k >= ival+1) then 1757 ebands%eig(ival+1:,ikpt,spin) = ebands%eig(ival+1:,ikpt,spin) + scissor_energy 1758 else 1759 write(msg,'(2a,4(a,i0))')& 1760 'Not enough bands to apply the scissor operator. ',ch10,& 1761 'spin: ',spin,' ikpt: ',ikpt,' nband_k: ',nband_k,' but valence index: ',ival 1762 ABI_COMMENT(msg) 1763 end if 1764 1765 end do 1766 end do 1767 1768 ! Recalculate the Fermi level and occupation factors. 1769 ! For Semiconductors only the Fermi level is changed (in the middle of the new gap) 1770 spinmagntarget_ = -99.99_dp !?; if (PRESENT(spinmagntarget)) spinmagntarget_=spinmagntarget 1771 call ebands_update_occ(ebands, spinmagntarget_) 1772 1773 end subroutine ebands_apply_scissors
m_ebands/ebands_calc_nelect [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_calc_nelect
FUNCTION
Compute nelect from Fermi level and Temperature.
INPUTS
OUTPUT
SOURCE
2747 real(dp) pure function ebands_calc_nelect(self, kt, fermie) result(nelect) 2748 2749 !Arguments ------------------------------------ 2750 !scalars 2751 class(ebands_t),intent(in) :: self 2752 real(dp),intent(in) :: kt, fermie 2753 2754 !Local variables------------------------------- 2755 !scalars 2756 integer :: spin, ik, ib 2757 real(dp) :: ofact 2758 2759 ! ************************************************************************* 2760 2761 ofact = two / (self%nsppol * self%nspinor) 2762 nelect = zero 2763 do spin=1,self%nsppol 2764 do ik=1,self%nkpt 2765 do ib=1,self%nband(ik + (spin-1)*self%nkpt) 2766 nelect = nelect + self%wtk(ik) * occ_fd(self%eig(ib,ik,spin), kt, fermie) 2767 end do 2768 end do 2769 end do 2770 2771 nelect = ofact * nelect 2772 2773 end function ebands_calc_nelect
m_ebands/ebands_chop [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_chop
FUNCTION
Return a new ebands_t object with a selected number of bands between bstart and bstop
INPUTS
SOURCE
4081 type(ebands_t) function ebands_chop(self, bstart, bstop) result(new) 4082 4083 !Arguments ------------------------------------ 4084 !scalars 4085 class(ebands_t),intent(in) :: self 4086 integer,intent(in) :: bstart, bstop 4087 4088 !Local variables ------------------------------ 4089 integer :: mband, nkpt, nsppol 4090 4091 ! ********************************************************************* 4092 4093 ABI_CHECK_IRANGE(bstart, 1, self%mband, "Invalid bstart") 4094 ABI_CHECK_IRANGE(bstop, 1, self%mband, "Invalid bstop") 4095 ABI_CHECK_ILEQ(bstart, bstop, "bstart should be <= bstop") 4096 4097 ! First copy the bands 4098 call ebands_copy(self, new) 4099 4100 ! Now chop them 4101 ABI_FREE(new%eig) 4102 ABI_FREE(new%occ) 4103 ABI_FREE(new%doccde) 4104 4105 mband = bstop - bstart + 1 4106 nkpt = self%nkpt 4107 nsppol = self%nsppol 4108 4109 ABI_MALLOC(new%eig, (mband, nkpt, nsppol)) 4110 ABI_MALLOC(new%occ, (mband, nkpt, nsppol)) 4111 ABI_MALLOC(new%doccde, (mband, nkpt, nsppol)) 4112 4113 new%mband = mband 4114 new%nband = mband 4115 new%eig = self%eig(bstart:bstop,:,:) 4116 new%occ = self%occ(bstart:bstop,:,:) 4117 new%doccde = self%doccde(bstart:bstop,:,:) 4118 4119 new%bantot = sum(new%nband) 4120 4121 end function ebands_chop
m_ebands/ebands_copy [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_copy
FUNCTION
This subroutine performs a deep copy of an ebands_t datatype.
INPUTS
ibands<ebands_t>=The data type to be copied.
OUTPUT
obands<ebands_t>=The copy.
SOURCE
1062 subroutine ebands_copy(ibands, obands) 1063 1064 !Arguments ------------------------------------ 1065 !scalars 1066 class(ebands_t),intent(in) :: ibands 1067 class(ebands_t),intent(out) :: obands 1068 1069 ! ********************************************************************* 1070 1071 call ebands_free(obands) 1072 1073 ! Copy scalars 1074 obands%bantot = ibands%bantot 1075 obands%ivalence = ibands%ivalence 1076 obands%mband = ibands%mband 1077 obands%nkpt = ibands%nkpt 1078 obands%nspinor = ibands%nspinor 1079 obands%nsppol = ibands%nsppol 1080 obands%occopt = ibands%occopt 1081 obands%kptopt = ibands%kptopt 1082 obands%nshiftk_orig = ibands%nshiftk_orig 1083 obands%nshiftk = ibands%nshiftk 1084 1085 obands%cellcharge = ibands%cellcharge 1086 obands%extrael = ibands%extrael 1087 obands%entropy = ibands%entropy 1088 obands%fermie = ibands%fermie 1089 obands%fermih = ibands%fermih 1090 obands%nelect = ibands%nelect 1091 obands%ne_qFD = ibands%ne_qFD 1092 obands%nh_qFD = ibands%nh_qFD 1093 obands%tphysel = ibands%tphysel 1094 obands%tsmear = ibands%tsmear 1095 1096 obands%kptrlatt_orig = ibands%kptrlatt_orig 1097 obands%kptrlatt = ibands%kptrlatt 1098 1099 ! Copy allocatable arrays 1100 ! integer 1101 call alloc_copy(ibands%istwfk, obands%istwfk) 1102 call alloc_copy(ibands%nband , obands%nband ) 1103 call alloc_copy(ibands%npwarr, obands%npwarr) 1104 1105 ! real 1106 call alloc_copy(ibands%kptns , obands%kptns ) 1107 call alloc_copy(ibands%eig , obands%eig ) 1108 call alloc_copy(ibands%occ , obands%occ ) 1109 call alloc_copy(ibands%doccde, obands%doccde) 1110 call alloc_copy(ibands%wtk , obands%wtk ) 1111 call alloc_copy(ibands%shiftk_orig, obands%shiftk_orig) 1112 call alloc_copy(ibands%shiftk, obands%shiftk) 1113 1114 if (allocated(ibands%linewidth)) call alloc_copy(ibands%linewidth, obands%linewidth) 1115 1116 end subroutine ebands_copy
m_ebands/ebands_downsample [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_downsample
FUNCTION
Return a new ebands_t object of type ebands_t with a coarser IBZ contained in the inititial one.
INPUTS
cryst<crystal_t>=Info on unit cell and symmetries. in_kptrlatt(3,3)=Defines the sampling of the "small" IBZ. Must be submesh of the "fine" mesh. in_nshiftk= Number of shifts in the coarse k-mesh in_shiftk(3, in_nshiftk) = Shifts of the coarse k-mesh
SOURCE
3956 type(ebands_t) function ebands_downsample(self, cryst, in_kptrlatt, in_nshiftk, in_shiftk) result(new) 3957 3958 !Arguments ------------------------------------ 3959 !scalars 3960 integer,intent(in) :: in_nshiftk 3961 type(ebands_t),intent(in) :: self 3962 type(crystal_t),intent(in) :: cryst 3963 !arrays 3964 integer,intent(in) :: in_kptrlatt(3,3) 3965 real(dp),intent(in) :: in_shiftk(3, in_nshiftk) 3966 3967 !Local variables------------------------------- 3968 !scalars 3969 integer,parameter :: sppoldbl1 = 1 3970 integer :: new_nkbz , timrev, bantot, new_nkibz, ik_ibz, ikf, spin,mband, comm 3971 real(dp) :: dksqmax 3972 character(len=500) :: msg 3973 !arrays 3974 integer,allocatable :: ibz_c2f(:,:) 3975 integer :: new_kptrlatt(3,3) 3976 integer,allocatable :: istwfk(:),nband(:,:),npwarr(:) 3977 real(dp),allocatable :: new_kbz(:,:), new_wtk(:), new_kibz(:,:), doccde(:), eig(:), occ(:) 3978 real(dp),allocatable :: doccde_3d(:,:,:), eig_3d(:,:,:), occ_3d(:,:,:), new_shiftk(:,:) 3979 3980 ! ********************************************************************* 3981 3982 comm = xmpi_comm_self 3983 3984 ! Find IBZ associated to the new mesh. 3985 call kpts_ibz_from_kptrlatt(cryst, in_kptrlatt, self%kptopt, in_nshiftk, in_shiftk, & 3986 new_nkibz, new_kibz, new_wtk, new_nkbz, new_kbz, new_kptrlatt=new_kptrlatt, new_shiftk=new_shiftk) 3987 3988 ! Costruct mapping IBZ_coarse --> IBZ_fine 3989 ! We don't change the value of nsppol hence sppoldbl1 is set to 1 3990 ABI_MALLOC(ibz_c2f, (new_nkibz*sppoldbl1, 6)) 3991 3992 timrev = kpts_timrev_from_kptopt(self%kptopt) 3993 call listkk(dksqmax, cryst%gmet, ibz_c2f, self%kptns, new_kibz, self%nkpt, new_nkibz, cryst%nsym, & 3994 sppoldbl1, cryst%symafm, cryst%symrel, timrev, comm, use_symrec=.False.) 3995 3996 if (dksqmax > tol12) then 3997 write(msg, '(a,es16.6,6a)' )& 3998 "At least one of the k-points could not be generated from a symmetrical one. dksqmax: ",dksqmax, ch10,& 3999 "kptrlatt of input ebands: ",trim(ltoa(pack(self%kptrlatt, mask=.True.))),ch10, & 4000 "downsampled K-mesh: ",trim(ltoa(pack(in_kptrlatt, mask=.True.))) 4001 ABI_ERROR(msg) 4002 end if 4003 4004 ABI_MALLOC(istwfk, (new_nkibz)) 4005 ABI_MALLOC(nband, (new_nkibz, self%nsppol)) 4006 ABI_MALLOC(npwarr, (new_nkibz)) 4007 4008 do ik_ibz=1,new_nkibz 4009 ikf = ibz_c2f(ik_ibz, 1) 4010 do spin=1,self%nsppol 4011 nband(ik_ibz, spin) = self%nband(ikf + (spin-1) * self%nkpt) 4012 end do 4013 istwfk(ik_ibz) = self%istwfk(ikf) 4014 npwarr(ik_ibz) = self%npwarr(ikf) 4015 end do 4016 4017 ! Recostruct eig, occ and doccde in the new IBZ. 4018 bantot = sum(nband); mband = maxval(nband) 4019 4020 ABI_MALLOC(doccde_3d, (mband, new_nkibz, self%nsppol)) 4021 ABI_MALLOC(eig_3d, (mband, new_nkibz, self%nsppol)) 4022 ABI_MALLOC(occ_3d, (mband, new_nkibz, self%nsppol)) 4023 4024 do spin=1,self%nsppol 4025 do ik_ibz=1,new_nkibz 4026 ikf = ibz_c2f(ik_ibz, 1) 4027 doccde_3d(:, ik_ibz, spin) = self%doccde(:, ikf, spin) 4028 eig_3d(:, ik_ibz, spin) = self%eig(:, ikf, spin) 4029 occ_3d(:, ik_ibz, spin) = self%occ(:, ikf, spin) 4030 end do 4031 end do 4032 4033 ! Have to pack data to call ebands_init (I wonder who decided to use vectors!) 4034 ABI_MALLOC(doccde, (bantot)) 4035 ABI_MALLOC(eig, (bantot)) 4036 ABI_MALLOC(occ, (bantot)) 4037 4038 call pack_eneocc(new_nkibz, self%nsppol, mband, nband, bantot, doccde_3d, doccde) 4039 call pack_eneocc(new_nkibz, self%nsppol, mband, nband, bantot, eig_3d, eig) 4040 call pack_eneocc(new_nkibz, self%nsppol, mband, nband, bantot, occ_3d, occ) 4041 4042 ABI_FREE(doccde_3d) 4043 ABI_FREE(eig_3d) 4044 ABI_FREE(occ_3d) 4045 4046 call ebands_init(bantot, new, self%nelect, self%ne_qFD, self%nh_qFD, self%ivalence, doccde, eig, istwfk, new_kibz, & 4047 nband, new_nkibz, npwarr, self%nsppol, self%nspinor, self%tphysel, self%tsmear, self%occopt, occ, new_wtk, & 4048 self%cellcharge, self%kptopt, in_kptrlatt, in_nshiftk, self%shiftk, new_kptrlatt, size(new_shiftk, dim=2), new_shiftk) 4049 4050 new%fermie = self%fermie 4051 new%fermih = self%fermih 4052 4053 ABI_FREE(istwfk) 4054 ABI_FREE(nband) 4055 ABI_FREE(npwarr) 4056 ABI_FREE(doccde) 4057 ABI_FREE(eig) 4058 ABI_FREE(occ) 4059 ABI_FREE(new_kibz) 4060 ABI_FREE(new_kbz) 4061 ABI_FREE(new_wtk) 4062 ABI_FREE(new_shiftk) 4063 ABI_FREE(ibz_c2f) 4064 4065 end function ebands_downsample
m_ebands/ebands_enclose_degbands [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_enclose_degbands
FUNCTION
Adjust ibmin and ibmax such that all the degenerate states are enclosed between ibmin and ibmax. The routine works for a given k-point a spin.
INPUTS
ebands<ebands_t>=The object describing the band structure. ikibz=Index of the k-point. spin=Spin index. tol_enedif=Tolerance on the energy difference.
OUTPUT
changed=.TRUE. if ibmin or ibmax has been changed.
[degblock(2,ndeg)]=Table allocated by the routine containing the index
of the bands in the `ndeg` degenerate sub-sets
degblock(1, ii) = first band index in the ii-th degenerate subset.
degblock(2, ii) = last band index in the ii-th degenerate subset.
SIDE EFFECTS
ibmin,ibmax=
Input: initial guess for the indices
Output: All the denerate states are between ibmin and ibmax
SOURCE
1866 subroutine ebands_enclose_degbands(ebands, ikibz, spin, ibmin, ibmax, changed, tol_enedif, degblock) 1867 1868 !Arguments ------------------------------------ 1869 !scalars 1870 class(ebands_t),intent(in) :: ebands 1871 integer,intent(in) :: ikibz,spin 1872 integer,intent(inout) :: ibmin,ibmax 1873 real(dp),intent(in) :: tol_enedif 1874 logical,intent(out) :: changed 1875 !arrays 1876 integer,allocatable,optional,intent(out) :: degblock(:,:) 1877 1878 !Local variables------------------------------- 1879 !scalars 1880 integer :: ib,ibmin_bkp,ibmax_bkp,ndeg 1881 real(dp) :: emin,emax 1882 1883 ! ************************************************************************* 1884 1885 ibmin_bkp = ibmin; ibmax_bkp = ibmax 1886 1887 emin = ebands%eig(ibmin,ikibz,spin) 1888 do ib=ibmin-1,1,-1 1889 if (ABS(ebands%eig(ib,ikibz,spin) - emin) > tol_enedif) then 1890 ibmin = ib +1 1891 EXIT 1892 else 1893 ibmin = ib 1894 end if 1895 end do 1896 1897 emax = ebands%eig(ibmax,ikibz,spin) 1898 do ib=ibmax+1,ebands%nband(ikibz+(spin-1)*ebands%nkpt) 1899 if ( ABS(ebands%eig(ib,ikibz,spin) - emax) > tol_enedif) then 1900 ibmax = ib - 1 1901 EXIT 1902 else 1903 ibmax = ib 1904 end if 1905 end do 1906 1907 changed = (ibmin /= ibmin_bkp) .or. (ibmax /= ibmax_bkp) 1908 1909 ! Compute degeneracy table. 1910 if (present(degblock)) then 1911 ! Count number of degeneracies. 1912 ndeg = 1 1913 do ib=ibmin+1,ibmax 1914 if ( abs(ebands%eig(ib,ikibz,spin) - ebands%eig(ib-1,ikibz,spin) ) > tol_enedif) ndeg = ndeg + 1 1915 end do 1916 ! Build degblock table. 1917 ABI_REMALLOC(degblock, (2, ndeg)) 1918 ndeg = 1; degblock(1, 1) = ibmin 1919 do ib=ibmin+1,ibmax 1920 if ( abs(ebands%eig(ib,ikibz,spin) - ebands%eig(ib-1,ikibz,spin) ) > tol_enedif) then 1921 degblock(2, ndeg) = ib - 1 1922 ndeg = ndeg + 1 1923 degblock(1, ndeg) = ib 1924 end if 1925 end do 1926 degblock(2, ndeg) = ibmax 1927 end if 1928 1929 end subroutine ebands_enclose_degbands
m_ebands/ebands_expandk [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_expandk
FUNCTION
Return a new object of type ebands_t corresponding to a list of k-points specified in input. Symmetry properties of the eigenvectors are used to symmetrize energies and occupation numbers.
INPUTS
inb<ebands_t>=Initial band structure with energies in the IBZ. ecut_eff=Effective cutoff energy i.e. ecut * dilatmx**2 force_istwfk1=If True, istwfk if forced to 1 for all the k-points in the BZ.
OUTPUT
dksqmax=maximal value of the norm**2 of the difference between
a kpt in the BZ and the closest k-point found in the inb%kpts set, using symmetries.
bz2ibz(nkpt2,6)=describe k point number of kpt1 that allows to
generate wavefunctions closest to given kpt2
bz2ibz(:,1)=k point number of kptns1
bz2ibz(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
(if 0, means no symmetry operation, equivalent to identity )
bz2ibz(:,3:5)=shift in reciprocal space to be given to kpt1a,
to give kpt1b, that is the closest to kpt2.
bz2ibz(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
outb<ebands_t>=band structure with energies in the BZ.
SOURCE
3801 subroutine ebands_expandk(inb, cryst, ecut_eff, force_istwfk1, dksqmax, bz2ibz, outb) 3802 3803 !Arguments ------------------------------------ 3804 !scalars 3805 real(dp),intent(in) :: ecut_eff 3806 real(dp),intent(out) :: dksqmax 3807 logical,intent(in) :: force_istwfk1 3808 class(ebands_t),intent(in) :: inb 3809 class(ebands_t),intent(out) :: outb 3810 class(crystal_t),intent(in) :: cryst 3811 !arrays 3812 integer,allocatable,intent(out) :: bz2ibz(:,:) 3813 3814 !Local variables------------------------------- 3815 !scalars 3816 integer,parameter :: istwfk_1=1,kptopt3=3 3817 integer :: nkfull,timrev,bantot,sppoldbl,npw_k,nsppol,istw 3818 integer :: ik_ibz,ikf,isym,itimrev,spin,mband,my_nkibz,comm 3819 logical :: isirred_k 3820 !character(len=500) :: msg 3821 !arrays 3822 integer :: g0(3) 3823 integer,allocatable :: istwfk(:),nband(:,:),npwarr(:),kg_k(:,:) 3824 real(dp),allocatable :: kfull(:,:),doccde(:),eig(:),occ(:),wtk(:),my_kibz(:,:) 3825 real(dp),allocatable :: doccde_3d(:,:,:),eig_3d(:,:,:),occ_3d(:,:,:) 3826 3827 ! ********************************************************************* 3828 3829 ABI_CHECK(inb%kptopt /= 0, "ebands_expandk does not support kptopt == 0") 3830 3831 comm = xmpi_comm_self 3832 nsppol = inb%nsppol 3833 3834 ! Note kptopt=3 3835 call kpts_ibz_from_kptrlatt(cryst, inb%kptrlatt, kptopt3, inb%nshiftk, inb%shiftk, & 3836 my_nkibz, my_kibz, wtk, nkfull, kfull) ! new_kptrlatt, new_shiftk) 3837 3838 ABI_FREE(my_kibz) 3839 ABI_FREE(wtk) 3840 3841 ! Costruct full BZ and create mapping BZ --> IBZ 3842 ! Note: 3843 ! - we don't change the value of nsppol hence sppoldbl is set to 1 3844 ! - we use symrel so that bz2ibz can be used to reconstruct the wavefunctions. 3845 ! 3846 sppoldbl = 1 !; if (any(cryst%symafm == -1) .and. inb%nsppol == 1) sppoldbl=2 3847 ABI_MALLOC(bz2ibz, (nkfull*sppoldbl,6)) 3848 3849 timrev = kpts_timrev_from_kptopt(inb%kptopt) 3850 call listkk(dksqmax, cryst%gmet, bz2ibz, inb%kptns, kfull, inb%nkpt, nkfull, cryst%nsym, & 3851 sppoldbl, cryst%symafm, cryst%symrel, timrev, comm, use_symrec=.False.) 3852 3853 ABI_MALLOC(wtk, (nkfull)) 3854 wtk = one / nkfull ! weights normalized to one 3855 3856 ABI_MALLOC(istwfk, (nkfull)) 3857 ABI_MALLOC(nband, (nkfull, nsppol)) 3858 ABI_MALLOC(npwarr, (nkfull)) 3859 3860 if (any(cryst%symrel(:,:,1) /= identity_3d) .and. any(abs(cryst%tnons(:,1)) > tol10) ) then 3861 ABI_ERROR('The first symmetry is not the identity operator!') 3862 end if 3863 3864 do ikf=1,nkfull 3865 ik_ibz = bz2ibz(ikf,1) 3866 isym = bz2ibz(ikf,2) 3867 itimrev = bz2ibz(ikf,6) 3868 g0 = bz2ibz(ikf,3:5) ! IS(k_ibz) + g0 = k_bz 3869 isirred_k = (isym == 1 .and. itimrev == 0 .and. all(g0 == 0)) 3870 3871 do spin=1,nsppol 3872 nband(ikf,spin) = inb%nband(ik_ibz+(spin-1)*inb%nkpt) 3873 end do 3874 3875 if (force_istwfk1) then 3876 call get_kg(kfull(:,ikf),istwfk_1,ecut_eff,cryst%gmet,npw_k,kg_k) 3877 ABI_FREE(kg_k) 3878 istwfk(ikf) = 1 3879 npwarr(ikf) = npw_k 3880 else 3881 if (isirred_k) then 3882 istwfk(ikf) = inb%istwfk(ik_ibz) 3883 npwarr(ikf) = inb%npwarr(ik_ibz) 3884 else 3885 istw = set_istwfk(kfull(:,ikf)) 3886 call get_kg(kfull(:,ikf),istw,ecut_eff,cryst%gmet,npw_k,kg_k) 3887 ABI_FREE(kg_k) 3888 istwfk(ikf) = istw 3889 npwarr(ikf) = npw_k 3890 end if 3891 end if 3892 end do 3893 3894 ! Recostruct eig, occ and doccde in the BZ. 3895 bantot = sum(nband); mband = maxval(nband) 3896 3897 ABI_MALLOC(doccde_3d, (mband, nkfull, nsppol)) 3898 ABI_MALLOC(eig_3d, (mband, nkfull, nsppol)) 3899 ABI_MALLOC(occ_3d, (mband, nkfull, nsppol)) 3900 3901 do spin=1,nsppol 3902 do ikf=1,nkfull 3903 ik_ibz = bz2ibz(ikf,1) 3904 doccde_3d(:,ikf,spin) = inb%doccde(:,ik_ibz,spin) 3905 eig_3d(:,ikf,spin) = inb%eig(:,ik_ibz,spin) 3906 occ_3d(:,ikf,spin) = inb%occ(:,ik_ibz,spin) 3907 end do 3908 end do 3909 3910 ! Have to pack data to call ebands_init (I wonder who decided to use vectors!) 3911 ABI_MALLOC(doccde, (bantot)) 3912 ABI_MALLOC(eig, (bantot)) 3913 ABI_MALLOC(occ, (bantot)) 3914 3915 call pack_eneocc(nkfull, nsppol, mband, nband, bantot, doccde_3d, doccde) 3916 call pack_eneocc(nkfull, nsppol, mband, nband, bantot, eig_3d, eig) 3917 call pack_eneocc(nkfull, nsppol, mband, nband, bantot, occ_3d, occ) 3918 3919 ABI_FREE(doccde_3d) 3920 ABI_FREE(eig_3d) 3921 ABI_FREE(occ_3d) 3922 3923 call ebands_init(bantot, outb, inb%nelect, inb%ne_qFD, inb%nh_qFD, inb%ivalence, doccde, eig, istwfk, kfull, & 3924 nband, nkfull, npwarr, nsppol, inb%nspinor, inb%tphysel, inb%tsmear, inb%occopt, occ, wtk, & 3925 inb%cellcharge, kptopt3, inb%kptrlatt_orig, inb%nshiftk_orig, inb%shiftk_orig, inb%kptrlatt, inb%nshiftk, inb%shiftk) 3926 3927 ABI_FREE(istwfk) 3928 ABI_FREE(nband) 3929 ABI_FREE(npwarr) 3930 ABI_FREE(doccde) 3931 ABI_FREE(eig) 3932 ABI_FREE(occ) 3933 ABI_FREE(wtk) 3934 ABI_FREE(kfull) 3935 3936 end subroutine ebands_expandk
m_ebands/ebands_free [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_free
FUNCTION
Deallocates the components of the ebands_t structured datatype
INPUTS
ebands<ebands_t>=The data type to be deallocated.
OUTPUT
Deallocate the dynamic arrays in the ebands_t type.
SOURCE
1022 subroutine ebands_free(ebands) 1023 1024 !Arguments ------------------------------------ 1025 class(ebands_t),intent(inout) :: ebands 1026 ! ************************************************************************* 1027 1028 ABI_SFREE(ebands%istwfk) 1029 ABI_SFREE(ebands%nband) 1030 ABI_SFREE(ebands%npwarr) 1031 ABI_SFREE(ebands%kptns) 1032 ABI_SFREE(ebands%eig) 1033 ABI_SFREE(ebands%linewidth) 1034 ABI_SFREE(ebands%occ) 1035 ABI_SFREE(ebands%doccde) 1036 ABI_SFREE(ebands%wtk) 1037 !ABI_SFREE(ebands%velocity) 1038 !ABI_SFREE(ebands%kTmesh) 1039 ABI_SFREE(ebands%shiftk_orig) 1040 ABI_SFREE(ebands%shiftk) 1041 1042 end subroutine ebands_free
m_ebands/ebands_from_dtset [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_from_dtset
FUNCTION
Build and return a new ebands_t datatype. Dimensions are taken from the abinit dataset.
INPUTS
dtset<dataset_type>=Abinit dataset npwarr(dtset%nkpt)=Number of G-vectors for each k-point. [nband]= If present, use these values instead of dtset%nband
OUTPUT
ebands<ebands_t>=The ebands_t datatype completely initialized. The Fermi level and the entropy are set to zero.
SOURCE
965 type(ebands_t) function ebands_from_dtset(dtset, npwarr, nband) result(new) 966 967 !Arguments ------------------------------------ 968 !scalars 969 type(dataset_type),target,intent(in) :: dtset 970 integer,target,optional,intent(in) :: nband(dtset%nkpt * dtset%nsppol) 971 !arrays 972 integer,intent(in) :: npwarr(dtset%nkpt) 973 974 !Local variables------------------------------- 975 !scalars 976 integer :: bantot 977 !arrays 978 real(dp),allocatable :: ugly_doccde(:), ugly_ene(:), ugly_occ(:) 979 integer,pointer :: nband__(:) 980 ! ************************************************************************* 981 982 nband__ => dtset%nband; if (present(nband)) nband__ => nband 983 984 ! Have to use ugly 1d vectors to call ebands_init 985 bantot = sum(nband__) 986 ABI_CALLOC(ugly_doccde, (bantot)) 987 ABI_CALLOC(ugly_ene, (bantot)) 988 ABI_CALLOC(ugly_occ, (bantot)) 989 990 call ebands_init(bantot, new, dtset%nelect, dtset%ne_qFD, dtset%nh_qFD, dtset%ivalence, ugly_doccde, ugly_ene, & 991 dtset%istwfk, dtset%kptns, nband__, dtset%nkpt, & 992 npwarr, dtset%nsppol, dtset%nspinor, dtset%tphysel, dtset%tsmear, dtset%occopt, ugly_occ, dtset%wtk,& 993 dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 994 dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 995 996 !new%extrael = dtset%eph_extrael 997 998 ABI_FREE(ugly_doccde) 999 ABI_FREE(ugly_ene) 1000 ABI_FREE(ugly_occ) 1001 1002 end function ebands_from_dtset
m_ebands/ebands_from_hdr [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_from_hdr
FUNCTION
This subroutine initializes the ebands_t datatype from the abinit header by calling the main creation method.
INPUTS
Hdr<hdr_type>=Abinit header. mband=Maximum number of bands. ene3d(mband,Hdr%nkpt,Hdr%nsppol)=Energies. [nelect]=Number of electrons per unit cell. Optional argument that can be used for performing a ridid shift of the fermi level. in the case of metallic occupancies. If not specified, nelect will be initialized from Hdr.
OUTPUT
ebands<ebands_t>=The ebands_t datatype completely initialized.
SOURCE
904 type(ebands_t) function ebands_from_hdr(hdr, mband, ene3d, nelect) result(ebands) 905 906 !Arguments ------------------------------------ 907 !scalars 908 integer,intent(in) :: mband 909 type(hdr_type),intent(in) :: hdr 910 real(dp),optional,intent(in) :: nelect 911 !arrays 912 real(dp),intent(in) :: ene3d(mband,hdr%nkpt,hdr%nsppol) 913 914 !Local variables------------------------------- 915 !scalars 916 real(dp) :: my_nelect 917 !arrays 918 real(dp),allocatable :: ugly_doccde(:),ugly_ene(:) 919 920 ! ************************************************************************* 921 922 my_nelect = hdr%nelect; if (present(nelect)) my_nelect = nelect 923 924 ! Have to use ugly 1d vectors to call ebands_init 925 ABI_CALLOC(ugly_doccde, (hdr%bantot)) 926 ABI_MALLOC(ugly_ene, (hdr%bantot)) 927 928 call pack_eneocc(hdr%nkpt, hdr%nsppol, mband, hdr%nband, hdr%bantot, ene3d, ugly_ene) 929 930 call ebands_init(hdr%bantot, ebands, my_nelect, hdr%ne_qFD, hdr%nh_qFD, hdr%ivalence, & 931 ugly_doccde, ugly_ene, hdr%istwfk, hdr%kptns, hdr%nband, hdr%nkpt, & 932 hdr%npwarr, hdr%nsppol, hdr%nspinor, hdr%tphysel, hdr%tsmear, hdr%occopt, hdr%occ, hdr%wtk, & 933 hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 934 935 ! Copy the fermi level reported in the header 936 ebands%fermie = hdr%fermie 937 ebands%fermih = hdr%fermih 938 939 ABI_FREE(ugly_doccde) 940 ABI_FREE(ugly_ene) 941 942 end function ebands_from_hdr
m_ebands/ebands_get_bandenergy [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_bandenergy
FUNCTION
Return the band energy (weighted sum of occupied eigenvalues)
INPUTS
OUTPUT
NOTES
TODO Likely this expression is not accurate since it is not variatonal One should use band_energy = \int e N(e) de for e<Ef , where N(e) is the e-DOS
SOURCE
1483 pure real(dp) function ebands_get_bandenergy(ebands) result(band_energy) 1484 1485 !Arguments ------------------------------------ 1486 !scalars 1487 class(ebands_t),intent(in) :: ebands 1488 1489 !Local variables------------------------------- 1490 integer :: spin,ikibz,nband_k 1491 real(dp) :: wtk 1492 ! ********************************************************************* 1493 1494 band_energy=zero 1495 do spin=1,ebands%nsppol 1496 do ikibz=1,ebands%nkpt 1497 wtk=ebands%wtk(ikibz) 1498 nband_k=ebands%nband(ikibz+(spin-1)*ebands%nkpt) 1499 band_energy = band_energy + wtk*SUM( ebands%eig(1:nband_k,ikibz,spin)*ebands%occ(1:nband_k,ikibz,spin) ) 1500 end do 1501 end do 1502 1503 end function ebands_get_bandenergy
m_ebands/ebands_get_bands_e0 [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_bands_e0
FUNCTION
Find min/max band indices crossing energy e0 min/max are returned in brange_spin(1:2, spin) for each spin. If no band crosses e0, bmin is set to +huge(1) and bmax to -huge(1) and ierr != 0
INPUTS
OUTPUT
SOURCE
1949 subroutine ebands_get_bands_e0(ebands, e0, brange_spin, ierr) 1950 1951 !Arguments ------------------------------------ 1952 !scalars 1953 class(ebands_t),intent(in) :: ebands 1954 real(dp),intent(in) :: e0 1955 integer,intent(out) :: brange_spin(2, ebands%nsppol) 1956 integer,intent(out) :: ierr 1957 1958 !Local variables------------------------------- 1959 integer :: band, spin, bmin, bmax 1960 real(dp) :: emin, emax 1961 1962 ! ************************************************************************* 1963 1964 ierr = 0 1965 do spin=1,ebands%nsppol 1966 bmin = +huge(1); bmax = -huge(1) 1967 1968 do band=1,minval(ebands%nband) 1969 emin = minval(ebands%eig(band, :, spin)) 1970 emax = maxval(ebands%eig(band, :, spin)) 1971 if (emin <= e0 .and. emax >= e0) then 1972 bmin = min(bmin, band) 1973 bmax = max(bmax, band) 1974 end if 1975 end do 1976 1977 brange_spin(:, spin) = [bmin, bmax] 1978 if (bmin == +huge(1)) ierr = ierr + 1 1979 end do 1980 1981 end subroutine ebands_get_bands_e0
m_ebands/ebands_get_bands_from_erange [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_bands_from_erange
FUNCTION
Return the indices of the min and max band index within an energy window.
INPUTS
elow, ehigh: Min and max energy
OUTPUT
bstart, bstop: Min and max band index. Initialized to bstart = huge(1); bstop = -huge(1)
SOURCE
1580 pure subroutine ebands_get_bands_from_erange(ebands, elow, ehigh, bstart, bstop) 1581 1582 !Arguments ------------------------------------ 1583 !scalars 1584 class(ebands_t),intent(in) :: ebands 1585 real(dp),intent(in) :: elow, ehigh 1586 integer,intent(out) :: bstart, bstop 1587 1588 !Local variables------------------------------- 1589 integer :: band, ik, spin 1590 1591 ! ************************************************************************* 1592 1593 bstart = huge(1); bstop = -huge(1) 1594 do spin=1,ebands%nsppol 1595 do ik=1,ebands%nkpt 1596 do band=1,ebands%nband(ik + (spin - 1) * ebands%nkpt) 1597 if (ebands%eig(band, ik , spin) >= elow .and. ebands%eig(band, ik , spin) <= ehigh) then 1598 bstart = min(bstart, band) 1599 bstop = max(bstop, band) 1600 end if 1601 end do 1602 end do 1603 end do 1604 1605 end subroutine ebands_get_bands_from_erange
m_ebands/ebands_get_carriers [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_carriers
FUNCTION
Compute number of electrons (e) and holes (h) per unit cell from a given list of `ntemp` temperatures `kTmesh` and chemical potentials `mu_e`. Return n_ehst(2, nsppol, ntemp) where the first dimension if for electrons/holes. If nsppol == 2, the second dimension is the number of e/h for spin else the total number of e/h summed over spins. To discern between electrons and holes in semiconductors we assume that ef is inside the gap.
SOURCE
6108 subroutine ebands_get_carriers(self, ntemp, kTmesh, mu_e, n_ehst) 6109 6110 !Arguments ------------------------------------ 6111 !scalars 6112 class(ebands_t),intent(in) :: self 6113 integer,intent(in) :: ntemp 6114 !arrays 6115 real(dp),intent(in) :: kTmesh(ntemp), mu_e(ntemp) 6116 real(dp),intent(out) :: n_ehst(2, self%nsppol, ntemp) 6117 6118 !Local variables------------------------------- 6119 integer :: spin, ik_ibz, ib, itemp 6120 real(dp) :: max_occ, wtk, eig_nk 6121 6122 ! ********************************************************************* 6123 6124 max_occ = two / (self%nspinor * self%nsppol) 6125 n_ehst = zero 6126 6127 do spin=1,self%nsppol 6128 do ik_ibz=1,self%nkpt 6129 wtk = self%wtk(ik_ibz) 6130 do ib=1, self%nband(ik_ibz + (spin-1) * self%nkpt) 6131 eig_nk = self%eig(ib, ik_ibz, spin) 6132 6133 do itemp=1,ntemp 6134 if (eig_nk >= mu_e(itemp)) then 6135 ! electron (assuming ef inside the gap if semiconductor) 6136 n_ehst(1, spin, itemp) = n_ehst(1, spin, itemp) + & 6137 wtk * occ_fd(eig_nk, kTmesh(itemp), mu_e(itemp)) * max_occ 6138 else 6139 ! holes 6140 n_ehst(2, spin, itemp) = n_ehst(2, spin, itemp) + & 6141 wtk * (one - occ_fd(eig_nk, kTmesh(itemp), mu_e(itemp))) * max_occ 6142 end if 6143 end do 6144 6145 end do 6146 end do 6147 end do 6148 6149 end subroutine ebands_get_carriers
m_ebands/ebands_get_edos [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_edos
FUNCTION
Calculate the electronic density of states from ebands_t
INPUTS
ebands<ebands_t>=Band structure object. cryst<cryst_t>=Info on the crystalline structure. intmeth= 1 for Gaussian, 2 or -2 for tetrahedra (-2 if Blochl corrections must be included). If nkpt == 1 (Gamma only), the routine fallbacks to gaussian method. step=Step on the linear mesh in Ha. If <0, the routine will use the mean of the energy level spacing broad=Gaussian broadening, If <0, the routine will use a default value for the broadening computed from the mean of the energy level spacing. No meaning for tetrahedra comm=MPI communicator
OUTPUT
edos<edos_t>=Electronic DOS and IDOS.
SOURCE
3146 type(edos_t) function ebands_get_edos(ebands, cryst, intmeth, step, broad, comm) result(edos) 3147 3148 !Arguments ------------------------------------ 3149 !scalars 3150 integer,intent(in) :: intmeth,comm 3151 real(dp),intent(in) :: step,broad 3152 class(ebands_t),target,intent(in) :: ebands 3153 type(crystal_t),intent(in) :: cryst 3154 3155 !Local variables------------------------------- 3156 !scalars 3157 integer :: nw,spin,band,ikpt,ief,ihf,nproc,my_rank,ierr,cnt,bcorr 3158 real(dp) :: max_ene,min_ene,wtk,max_occ 3159 character(len=500) :: msg 3160 type(htetra_t) :: tetra 3161 !arrays 3162 real(dp) :: eminmax_spin(2,ebands%nsppol) 3163 real(dp),allocatable :: wme0(:),wdt(:,:),tmp_eigen(:) 3164 3165 ! ********************************************************************* 3166 3167 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 3168 ierr = 0 3169 3170 edos%nkibz = ebands%nkpt; edos%nsppol = ebands%nsppol; edos%nspinor = ebands%nspinor 3171 edos%intmeth = intmeth 3172 edos%nelect = ebands%nelect 3173 3174 if (ebands%nkpt == 1) then 3175 ABI_COMMENT("Cannot use tetrahedra for e-DOS when nkpt == 1. Switching to gaussian method") 3176 edos%intmeth = 1 3177 end if 3178 3179 edos%broad = broad; edos%step = step 3180 3181 ! Compute the linear mesh so that it encloses all bands. 3182 eminmax_spin = ebands_get_minmax(ebands, "eig") 3183 min_ene = minval(eminmax_spin(1, :)); min_ene = min_ene - 0.1_dp * abs(min_ene) 3184 max_ene = maxval(eminmax_spin(2, :)); max_ene = max_ene + 0.1_dp * abs(max_ene) 3185 3186 nw = nint((max_ene - min_ene) / edos%step) + 1; edos%nw = nw 3187 3188 ABI_MALLOC(edos%mesh, (nw)) 3189 edos%mesh = arth(min_ene, edos%step, nw) 3190 3191 ABI_CALLOC(edos%gef, (0:edos%nsppol)) 3192 ABI_CALLOC(edos%ghf, (0:edos%nsppol)) 3193 ABI_CALLOC(edos%dos, (nw, 0:edos%nsppol)) 3194 ABI_CALLOC(edos%idos, (nw, 0:edos%nsppol)) 3195 3196 select case (edos%intmeth) 3197 case (1) 3198 !call wrtout(std_out, " Computing electron-DOS with Gaussian method") 3199 !call wrtout(std_out, sjoin(" broadening: ", ftoa(edos%broad * Ha_eV), " (eV), step: ", ftoa(edos%step * Ha_eV), "(eV), npts: ", 3200 !itoa(nw))) 3201 ! Gaussian 3202 ABI_MALLOC(wme0, (nw)) 3203 cnt = 0 3204 do spin=1,edos%nsppol 3205 do ikpt=1,ebands%nkpt 3206 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism 3207 wtk = ebands%wtk(ikpt) 3208 do band=1,ebands%nband(ikpt+(spin-1)*ebands%nkpt) 3209 wme0 = edos%mesh - ebands%eig(band, ikpt, spin) 3210 edos%dos(:, spin) = edos%dos(:, spin) + wtk * gaussian(wme0, edos%broad) 3211 end do 3212 end do 3213 end do 3214 ABI_FREE(wme0) 3215 call xmpi_sum(edos%dos, comm, ierr) 3216 3217 case (2, -2) 3218 !call wrtout(std_out, " Computing electron-DOS with tetrahedron method") 3219 ! Consistency test 3220 if (any(ebands%nband /= ebands%nband(1)) ) ABI_ERROR('for tetrahedra, nband(:) must be constant') 3221 3222 ! Build tetra object. 3223 tetra = tetra_from_kptrlatt(cryst, ebands%kptopt, ebands%kptrlatt, & 3224 ebands%nshiftk, ebands%shiftk, ebands%nkpt, ebands%kptns, comm, msg, ierr) 3225 ABI_CHECK(ierr == 0, msg) 3226 3227 ! For each spin and band, interpolate over kpoints, calculate integration weights and DOS contribution. 3228 ABI_MALLOC(tmp_eigen, (ebands%nkpt)) 3229 ABI_MALLOC(wdt, (nw, 2)) 3230 3231 bcorr = 0; if (intmeth == -2) bcorr = 1 3232 cnt = 0 3233 do spin=1,ebands%nsppol 3234 do band=1,ebands%nband(1) 3235 ! For each band get its contribution 3236 tmp_eigen = ebands%eig(band,:,spin) 3237 do ikpt=1,ebands%nkpt 3238 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 3239 3240 ! Calculate integration weights at each irred k-point (Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]) 3241 call tetra%get_onewk(ikpt, bcorr, nw, ebands%nkpt, tmp_eigen, min_ene, max_ene, one, wdt) 3242 3243 edos%dos(:,spin) = edos%dos(:,spin) + wdt(:, 1) * ebands%wtk(ikpt) 3244 edos%idos(:,spin) = edos%idos(:,spin) + wdt(:, 2) * ebands%wtk(ikpt) 3245 end do ! ikpt 3246 end do ! band 3247 end do ! spin 3248 3249 call xmpi_sum(edos%dos, comm, ierr) 3250 call xmpi_sum(edos%idos, comm, ierr) 3251 3252 ! Free memory 3253 ABI_FREE(tmp_eigen) 3254 ABI_FREE(wdt) 3255 call tetra%free() 3256 3257 ! Filter so that dos[i] is always >= 0 and idos is monotonic 3258 ! IDOS is computed afterwards with simpson 3259 where (edos%dos(:,1:) <= zero) 3260 edos%dos(:,1:) = zero 3261 end where 3262 3263 case default 3264 ABI_ERROR(sjoin("Wrong integration method:", itoa(intmeth))) 3265 end select 3266 3267 ! Compute total DOS and IDOS by summing the two spin channels. 3268 max_occ = two / (ebands%nspinor * ebands%nsppol) 3269 edos%dos(:, 0) = max_occ * sum(edos%dos(:,1:), dim=2) 3270 3271 if (edos%intmeth == 1) then 3272 do spin=1,edos%nsppol 3273 call simpson_int(nw, edos%step, edos%dos(:,spin), edos%idos(:,spin)) 3274 end do 3275 end if 3276 edos%idos(:, 0) = max_occ * sum(edos%idos(:,1:), dim=2) 3277 3278 ! Use bisection to find the Fermi level at T = 0 3279 ! Warning: this code assumes idos[i+1] >= idos[i]. This condition may not be 3280 ! fullfilled if we use tetra and this is the reason why we have filtered the DOS. 3281 if (ebands%occopt == 9) then 3282 ihf = bisect(edos%idos(:,0), ebands%nelect-ebands%nh_qFD) 3283 ief = bisect(edos%idos(:,0), ebands%nelect+ebands%ne_qFD) 3284 else 3285 ief = bisect(edos%idos(:,0), ebands%nelect) 3286 ihf = ief 3287 end if 3288 3289 ! Handle out of range condition. 3290 if (ief == 0 .or. ief == nw) then 3291 write(msg,"(3a)")& 3292 "Bisection could not find an initial guess for the Fermi level!",ch10,& 3293 "Possible reasons: not enough bands or wrong number of electrons" 3294 ABI_WARNING(msg) 3295 return 3296 else if (ihf == 0 .or. ihf == nw) then 3297 write(msg,"(3a)")& 3298 "Bisection could not find an initial guess for the holes Fermi level!",ch10,& 3299 "Possible reasons: not enough bands or wrong number of holes" 3300 ABI_WARNING(msg) 3301 return 3302 end if 3303 3304 ! TODO: Use linear interpolation to find an improved estimate of the Fermi level? 3305 edos%ief = ief 3306 edos%ihf = ihf 3307 do spin=0,edos%nsppol 3308 edos%gef(spin) = edos%dos(ief,spin) 3309 edos%ghf(spin) = edos%dos(ihf,spin) 3310 end do 3311 3312 !write(std_out,*)"fermie from ebands: ",ebands%fermie 3313 !write(std_out,*)"fermie from IDOS: ",edos%mesh(ief) 3314 !write(std_out,*)"gef:from ebands%fermie: " ,edos%dos(bisect(edos%mesh, ebands%fermie), 0) 3315 !write(std_out,*)"gef:from edos: " ,edos%gef(0) 3316 3317 end function ebands_get_edos
m_ebands/ebands_get_edos_matrix_elements [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_edos_matrix_elements
FUNCTION
Compute electron DOS and weighted e-DOS with weights given by precomputed scalar, vectorial or tensorial matrix elements. Weights are provided in input as (..., num_entries, bsize, nkpt, nsppol) tables, see below.
INPUTS
ebands<ebands_t>=Band structure object. cryst<cryst_t>=Info on the crystalline structure. bsize=Number of bands in bks_vals, bks_vecs and bks_tens Not necessarily equal to ebands%mband when brange is used. nvals=Number of scalar entries. Maybe zero bks_vals=Scalar matrix elements nvecs=Number of 3d-vectorial entries. Maybe zero bks_vecs=Vectorial matrix elements in Cartesian Coordinates ntens=Numer of 3x3 tensorial entries in Cartesian coordinates. Maybe zero bks_tens= Tensorial matrix elements (3x3) in Cartesian Coordinates intmeth= 1 for Gaussian, 2 or -2 for tetrahedra (-2 if Blochl corrections must be included). If nkpt == 1 (Gamma only), the routine fallbacks to the Gaussian method. step=Step on the linear mesh in Ha. If < 0, the routine will use the mean of the energy level spacing broad=Gaussian broadening, If < 0, the routine will use a default value for the broadening computed from the mean of the energy level spacing. No meaning for tetrahedra comm=MPI communicator [brange(2)]=Minimum and maximum band index. Default if not present is `full band range`. If given bsize must be equal: to brange(2) - brange(1) + 1 [erange(2)]=Minimum and maximum energy to be considered. Default if not present is `full energy range`.
OUTPUT
out_valsdos: (nw, 2, nvals, nsppol) array with DOS for scalar quantities if nvals > 0 out_vecsdos: (nw, 2, 3, nvecs, nsppol)) array with DOS weighted by vectorial terms if nvecs > 0 out_tensdos: (nw, 2,3, 3, ntens, nsppol) array with DOS weighted by tensorial terms if ntens > 0 All these arrays are allocated by the routine. The number of points is available in edos%nw. (nw, 1, ...) stores the weighted DOS (w-DOS) (nw, 2, ...) stores the integrated w-DOS
SOURCE
4527 type(edos_t) function ebands_get_edos_matrix_elements(ebands, cryst, bsize, & 4528 nvals, bks_vals, nvecs, bks_vecs, ntens, bks_tens, & 4529 intmeth, step, broad, out_valsdos, out_vecsdos, out_tensdos, comm, & 4530 brange, erange) result(edos) 4531 4532 !Arguments ------------------------------------ 4533 !scalars 4534 integer,intent(in) :: bsize, nvals, nvecs, ntens, intmeth, comm 4535 real(dp),intent(in) :: step, broad 4536 class(ebands_t),intent(in) :: ebands 4537 type(crystal_t),intent(in) :: cryst 4538 !arrays 4539 integer,optional,intent(in) :: brange(2) 4540 real(dp),optional,intent(in) :: erange(2) 4541 real(dp),intent(in) :: bks_vals(nvals, bsize, ebands%nkpt, ebands%nsppol) 4542 real(dp),intent(in) :: bks_vecs(3, nvecs, bsize, ebands%nkpt, ebands%nsppol) 4543 real(dp),intent(in) :: bks_tens(3, 3, ntens, bsize, ebands%nkpt, ebands%nsppol) 4544 real(dp),allocatable,intent(out) :: out_valsdos(:,:,:,:), out_vecsdos(:,:,:,:,:), out_tensdos(:,:,:,:,:,:) 4545 4546 !Local variables------------------------------- 4547 !scalars 4548 integer :: nproc, my_rank, nw, spin, band, ib, ik_ibz, cnt, idat, ierr, bcorr, time_opt 4549 integer :: ii, jj, ief, ihf, bmin_, bmax_ 4550 real(dp),parameter :: max_occ1 = one 4551 real(dp) :: emax, emin, wtk, max_occ 4552 real(dp) :: cpu, wall, gflops 4553 logical :: check_erange 4554 character(len=500) :: msg 4555 type(htetra_t) :: tetra 4556 !arrays 4557 real(dp) :: eminmax_spin(2,ebands%nsppol), vsum(3), tsum(3,3) 4558 real(dp),allocatable :: wme0(:),tmp_eigen(:), weights(:,:) 4559 4560 ! ********************************************************************* 4561 4562 call cwtime(cpu, wall, gflops, "start") 4563 4564 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 4565 4566 edos%nkibz = ebands%nkpt; edos%nsppol = ebands%nsppol; edos%nspinor = ebands%nspinor 4567 edos%intmeth = intmeth 4568 edos%nelect = ebands%nelect 4569 if (ebands%nkpt == 1) then 4570 ABI_COMMENT("Cannot use tetrahedra for e-DOS when nkpt == 1. Switching to gaussian method") 4571 edos%intmeth = 1 4572 end if 4573 4574 edos%broad = broad; edos%step = step 4575 4576 ! Define band range. 4577 bmin_ = 1; bmax_ = ebands%mband 4578 if (present(brange)) then 4579 bmin_ = brange(1); bmax_ = brange(2) 4580 end if 4581 4582 ABI_CHECK_IRANGE(bmin_, 1, ebands%mband, "Wrong bmin:") 4583 ABI_CHECK_IRANGE(bmax_, bmin_, ebands%mband, "Wrong bmax:") 4584 ABI_CHECK_IRANGE(bsize, 1, ebands%mband, "Wrong bsize:") 4585 4586 if (present(erange)) then 4587 ! use optional args if provided. 4588 emin = erange(1) 4589 emax = erange(2) 4590 check_erange = .True. 4591 else 4592 ! Compute the linear mesh so that it encloses all bands. 4593 eminmax_spin = ebands_get_minmax(ebands, "eig") 4594 emin = minval(eminmax_spin(1, :)); emin = emin - 0.1_dp * abs(emin) 4595 emax = maxval(eminmax_spin(2, :)); emax = emax + 0.1_dp * abs(emax) 4596 check_erange = .False. 4597 end if 4598 4599 nw = nint((emax - emin) / edos%step) + 1 4600 edos%nw = nw 4601 4602 ABI_MALLOC(edos%mesh, (nw)) 4603 edos%mesh = arth(emin, edos%step, nw) 4604 4605 ABI_CALLOC(edos%gef, (0:edos%nsppol)) 4606 ABI_CALLOC(edos%ghf, (0:edos%nsppol)) 4607 ABI_CALLOC(edos%dos, (nw, 0:edos%nsppol)) 4608 ABI_CALLOC(edos%idos, (nw, 0:edos%nsppol)) 4609 4610 ! Allocate output arrays depending on input. 4611 if (nvals > 0) then 4612 ABI_CALLOC(out_valsdos, (nw, 2, nvals, ebands%nsppol)) 4613 endif 4614 if (nvecs > 0) then 4615 ABI_CALLOC(out_vecsdos, (nw, 2, 3, nvecs, ebands%nsppol)) 4616 end if 4617 if (ntens > 0) then 4618 ABI_CALLOC(out_tensdos, (nw, 2, 3, 3, ntens, ebands%nsppol)) 4619 end if 4620 4621 time_opt = 0 ! This to preserve the previous behaviour in which TR was not used. 4622 4623 !call wrtout(std_out, " Computing DOS weighted by matrix elements.") 4624 select case (intmeth) 4625 case (1) 4626 ! Gaussian 4627 ABI_MALLOC(wme0, (nw)) 4628 cnt = 0 4629 do spin=1,ebands%nsppol 4630 do ik_ibz=1,ebands%nkpt 4631 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism 4632 wtk = ebands%wtk(ik_ibz) 4633 do band=bmin_,bmax_ 4634 ib = band - bmin_ + 1 4635 4636 if (check_erange) then 4637 if (ebands%eig(band, ik_ibz, spin) < emin - five * broad) cycle 4638 if (ebands%eig(band, ik_ibz, spin) > emax + five * broad) cycle 4639 end if 4640 4641 wme0 = edos%mesh - ebands%eig(band, ik_ibz, spin) 4642 wme0 = gaussian(wme0, broad) * wtk 4643 edos%dos(:,spin) = edos%dos(:,spin) + wme0(:) 4644 4645 ! scalars 4646 do idat=1,nvals 4647 out_valsdos(:, 1, idat, spin) = out_valsdos(:,1, idat, spin) + wme0(:) * bks_vals(idat, ib, ik_ibz, spin) 4648 ! FIXME: This is quite inefficient! Integration should be performed outside! 4649 call simpson_int(nw, step, out_valsdos(:,1, idat, spin), out_valsdos(:,2,idat,spin)) 4650 end do 4651 4652 ! vectors 4653 do idat=1,nvecs 4654 ! get components, symmetrize and accumulate. 4655 vsum = cryst%symmetrize_cart_vec3(bks_vecs(:, idat, ib, ik_ibz, spin), time_opt) 4656 do ii=1,3 4657 out_vecsdos(:, 1, ii, idat, spin) = out_vecsdos(:, 1, ii, idat, spin) + wme0(:) * vsum(ii) 4658 call simpson_int(nw, step, out_vecsdos(:,1,ii,idat,spin), out_vecsdos(:,2,ii,idat,spin)) 4659 end do 4660 end do 4661 4662 ! tensor 4663 do idat=1,ntens 4664 ! get components, symmetrize and accumulate. 4665 tsum = cryst%symmetrize_cart_tens33(bks_tens(:, :, idat, ib, ik_ibz, spin), time_opt) 4666 do ii=1,3 4667 do jj=1,3 4668 out_tensdos(:,1,jj,ii,idat,spin) = out_tensdos(:,1,jj,ii,idat,spin) + wme0(:) * tsum(jj,ii) 4669 call simpson_int(nw, step, out_tensdos(:,1,jj,ii,idat,spin), out_tensdos(:,2,jj,ii,idat,spin)) 4670 end do 4671 end do 4672 end do 4673 4674 end do !band 4675 end do !ik_ibz 4676 end do !spin 4677 4678 ABI_FREE(wme0) 4679 call xmpi_sum(edos%dos, comm, ierr) 4680 if (nvals > 0) call xmpi_sum(out_valsdos, comm, ierr) 4681 if (nvecs > 0) call xmpi_sum(out_vecsdos, comm, ierr) 4682 if (ntens > 0) call xmpi_sum(out_tensdos, comm, ierr) 4683 4684 case (2, -2) 4685 ! Consistency test 4686 ABI_CHECK(all(ebands%nband == ebands%nband(1)), 'For tetrahedra, nband(:) must be constant') 4687 4688 ! Build tetra object. 4689 tetra = tetra_from_kptrlatt(cryst, ebands%kptopt, ebands%kptrlatt, & 4690 ebands%nshiftk, ebands%shiftk, ebands%nkpt, ebands%kptns, comm, msg, ierr) 4691 ABI_CHECK(ierr == 0, msg) 4692 4693 ! For each spin and band, interpolate over kpoints, 4694 ! calculate integration weights and DOS contribution. 4695 ABI_MALLOC(tmp_eigen, (ebands%nkpt)) 4696 ABI_MALLOC(weights, (nw, 2)) 4697 4698 ! Blochl's corrections? 4699 bcorr = 0; if (intmeth == -2) bcorr = 1 4700 4701 cnt = 0 4702 do spin=1,ebands%nsppol 4703 do band=bmin_,bmax_ 4704 ! For each band get its contribution 4705 tmp_eigen = ebands%eig(band,:,spin) 4706 ib = band - bmin_ + 1 4707 4708 if (check_erange) then 4709 if (all(tmp_eigen < emin)) cycle 4710 if (all(tmp_eigen > emax)) cycle 4711 end if 4712 4713 do ik_ibz=1,ebands%nkpt 4714 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism 4715 4716 call tetra%get_onewk(ik_ibz, bcorr, nw, ebands%nkpt, tmp_eigen, emin, emax, max_occ1, weights) 4717 weights = weights * ebands%wtk(ik_ibz) 4718 4719 ! Compute DOS and IDOS 4720 edos%dos(:,spin) = edos%dos(:,spin) + weights(:, 1) 4721 edos%idos(:,spin) = edos%idos(:,spin) + weights(:, 2) 4722 4723 ! scalar 4724 !$OMP PARALLEL DO 4725 do idat=1,nvals 4726 out_valsdos(:, :, idat, spin) = out_valsdos(:, :, idat, spin) + weights(:, :) * bks_vals(idat, ib, ik_ibz, spin) 4727 end do 4728 4729 ! vector 4730 !$OMP PARALLEL DO PRIVATE(vsum) 4731 do idat=1,nvecs 4732 ! get components, symmetrize and accumulate. 4733 vsum = cryst%symmetrize_cart_vec3(bks_vecs(:, idat, ib, ik_ibz, spin), time_opt) 4734 do ii=1,3 4735 out_vecsdos(:, :, ii, idat, spin) = out_vecsdos(:, :, ii, idat, spin) + weights(:, :) * vsum(ii) 4736 end do 4737 end do 4738 4739 ! tensor 4740 !$OMP PARALLEL DO PRIVATE(tsum) 4741 do idat=1,ntens 4742 ! get components, symmetrize and accumulate. 4743 tsum = cryst%symmetrize_cart_tens33(bks_tens(:, :, idat, ib, ik_ibz, spin), time_opt) 4744 do ii=1,3 4745 do jj=1,3 4746 out_tensdos(:, :, jj, ii, idat, spin) = out_tensdos(:, :, jj, ii, idat, spin) + weights(:, :) * tsum(jj,ii) 4747 end do 4748 end do 4749 end do 4750 4751 end do ! ik_ibz 4752 end do ! band 4753 end do ! spin 4754 4755 ! Free memory 4756 ABI_FREE(weights) 4757 ABI_FREE(tmp_eigen) 4758 call tetra%free() 4759 4760 call xmpi_sum(edos%dos, comm, ierr) 4761 call xmpi_sum(edos%idos, comm, ierr) 4762 if (nvals > 0) call xmpi_sum(out_valsdos, comm, ierr) 4763 if (nvecs > 0) call xmpi_sum(out_vecsdos, comm, ierr) 4764 if (ntens > 0) call xmpi_sum(out_tensdos, comm, ierr) 4765 4766 case default 4767 ABI_ERROR(sjoin("Wrong integration method:", itoa(intmeth))) 4768 end select 4769 4770 ! Compute total DOS and IDOS 4771 max_occ = two / (ebands%nspinor * ebands%nsppol) 4772 edos%dos(:, 0) = max_occ * sum(edos%dos(:,1:), dim=2) 4773 4774 if (intmeth == 1) then 4775 do spin=1,edos%nsppol 4776 call simpson_int(nw, edos%step, edos%dos(:,spin), edos%idos(:,spin)) 4777 end do 4778 end if 4779 edos%idos(:, 0) = max_occ * sum(edos%idos(:,1:), dim=2) 4780 4781 ! Use bisection to find the Fermi level. 4782 ! Warning: this code assumes idos[i+1] >= idos[i]. This condition may not be 4783 ! fullfilled if we use tetra and this is the reason why we have filtered the DOS. 4784 if (ebands%occopt == 9) then 4785 ihf = bisect(edos%idos(:,0), ebands%nelect-ebands%nh_qFD) 4786 ief = bisect(edos%idos(:,0), ebands%nelect+ebands%ne_qFD) 4787 else 4788 ief = bisect(edos%idos(:,0), ebands%nelect) 4789 ihf = ief 4790 end if 4791 4792 ! Handle out of range condition. 4793 if (ief == 0 .or. ief == nw) then 4794 write(msg,"(a, f14.2, 4a)") & 4795 "Bisection could not find an initial guess for the Fermi level with nelect: ",ebands%nelect, ch10, & 4796 "Possible reasons: not enough bands for DOS or wrong number of electrons.", ch10, & 4797 "Returning from ebands_get_edos_matrix_elements without setting edos%ief !" 4798 ABI_WARNING(msg) 4799 return 4800 end if 4801 4802 if (ihf == 0 .or. ihf == nw) then 4803 write(msg,"(3a)")& 4804 "Bisection could not find an initial guess for the holes Fermi level!",ch10,& 4805 "Possible reasons: not enough bands or wrong number of holes" 4806 ABI_WARNING(msg) 4807 return 4808 end if 4809 4810 ! TODO: Use linear interpolation to find an improved estimate of the Fermi level? 4811 edos%ief = ief 4812 edos%ihf = ihf 4813 do spin=0,edos%nsppol 4814 edos%gef(spin) = edos%dos(ief,spin) 4815 edos%ghf(spin) = edos%dos(ihf,spin) 4816 end do 4817 4818 call cwtime_report(" ebands_get_edos_matrix_elements", cpu, wall, gflops) 4819 4820 end function ebands_get_edos_matrix_elements
m_ebands/ebands_get_erange [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_erange
FUNCTION
Compute the minimum and maximum energy enclosing a list of states specified by k-points and band indices.
INPUTS
ebands<ebands_t>=The object describing the band structure. nkpts=Number of k-points kpoints(3,nkpts)=K-points band_block(2,nkpts)=Gives for each k-points, the initial and the final band index to include.
OUTPUT
emin,emax=min and max energy
SOURCE
2005 subroutine ebands_get_erange(ebands, nkpts, kpoints, band_block, emin, emax) 2006 2007 !Arguments ------------------------------------ 2008 !scalars 2009 integer,intent(in) :: nkpts 2010 real(dp),intent(out) :: emin,emax 2011 class(ebands_t),intent(in) :: ebands 2012 !arrays 2013 integer,intent(in) :: band_block(2,nkpts) 2014 real(dp),intent(in) :: kpoints(3,nkpts) 2015 2016 !Local variables------------------------------- 2017 !scalars 2018 integer :: spin,ik,ikpt,cnt 2019 type(krank_t) :: krank 2020 2021 ! ************************************************************************* 2022 2023 krank = krank_new(ebands%nkpt, ebands%kptns) 2024 emin = huge(one); emax = -huge(one); cnt = 0 2025 2026 do spin=1,ebands%nsppol 2027 do ik=1,nkpts 2028 ikpt = krank%get_index(kpoints(:,ik)) 2029 if (ikpt == -1) then 2030 ABI_WARNING(sjoin("Cannot find k-point:", ktoa(kpoints(:,ik)))) 2031 cycle 2032 end if 2033 if (.not. (band_block(1,ik) >= 1 .and. band_block(2,ik) <= ebands%mband)) cycle 2034 cnt = cnt + 1 2035 emin = min(emin, minval(ebands%eig(band_block(1,ik):band_block(2,ik), ikpt, spin))) 2036 emax = max(emax, maxval(ebands%eig(band_block(1,ik):band_block(2,ik), ikpt, spin))) 2037 end do 2038 end do 2039 2040 call krank%free() 2041 2042 ! This can happen if wrong input. 2043 if (cnt == 0) then 2044 ABI_WARNING("None of the k-points/bands provided was found in ebands%") 2045 emin = minval(ebands%eig); emax = maxval(ebands%eig) 2046 end if 2047 2048 end subroutine ebands_get_erange
m_ebands/ebands_get_gaps [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_gaps
FUNCTION
Returns a gaps_t object with info on the fundamental and direct gap.
INPUTS
ebands<ebands_t>=Info on the band structure, the smearing technique and the physical temperature used.
OUTPUT
ierr=Return code (!=0 signals failure) gaps<gaps_t>=object with info on the gaps (caller is responsible for freeing the object).
SOURCE
436 type(gaps_t) function ebands_get_gaps(ebands, ierr) result(gaps) 437 438 !Arguments ------------------------------------ 439 !scalars 440 class(ebands_t),target,intent(in) :: ebands 441 integer,intent(out) :: ierr 442 443 !Local variables------------------------------- 444 !scalars 445 integer,parameter :: occopt3 = 3, prtvol0 = 0 446 real(dp),parameter :: spinmagntarget_ = -99.99_dp 447 real(dp) :: tsmear 448 type(ebands_t) :: tmp_ebands 449 !character(len=500) :: msg 450 451 ! ********************************************************************* 452 453 gaps = get_gaps_(ebands, ierr) 454 455 if (ierr /= 0) then 456 ! get_gaps_ will fail if we have a real metal/semimetal 457 ! but it's also possible to have a false negative if the input ebands represents a: 458 ! 459 ! 1) highly degenerate doped semiconductor with the Fermi level in the bands. 460 ! 2) Small gap semiconductor at relatively high T. 461 ! 462 ! Here I try to compute the gaps of an intrinsic semiconductor at low T with Fermi-Dirac. 463 ! This might still fail though and the caller should handle that. 464 call gaps%free() 465 call ebands_copy(ebands, tmp_ebands) 466 tsmear = 0.01_dp * eV_Ha 467 call ebands_set_scheme(tmp_ebands, occopt3, tsmear, spinmagntarget_, prtvol0, update_occ=.False.) 468 ! Remove extrael to go back to intrinsic system 469 if (ebands%extrael /= zero) tmp_ebands%nelect = ebands%nelect - ebands%extrael 470 !if (ebands%cellcharge /= zero) tmp_ebands%nelect = ebands%nelect + ebands%cellcharge 471 call ebands_update_occ(tmp_ebands, spinmagntarget_) 472 473 ! Try to compute gaps the again with new Fermi level at FD T = tsmear computed from update_occ. 474 ! Return ierr 475 gaps = get_gaps_(tmp_ebands, ierr) 476 call ebands_free(tmp_ebands) 477 end if 478 479 end function ebands_get_gaps
m_ebands/ebands_get_jdos [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_jdos
FUNCTION
Compute the joint density of states.
INPUTS
ebands<ebands_t>=Band structure object. cryst<cryst_t>=Info on the crystalline structure. intmeth= 1 for gaussian, 2 or -2 for tetrahedra (-2 if Blochl corrections must be included). step=Step on the linear mesh in Ha. If <0, the routine will use the mean of the energy level spacing broad=Gaussian broadening, If <0, the routine will use a default value for the broadening computed from the mean of the energy level spacing. No meaning if tetra method comm=MPI communicator
OUTPUT
SOURCE
4846 type(jdos_t) function ebands_get_jdos(ebands, cryst, intmeth, step, broad, comm, ierr) result (jdos) 4847 4848 !Arguments ------------------------------------ 4849 !scalars 4850 integer,intent(in) :: intmeth,comm 4851 integer,intent(out) :: ierr 4852 real(dp),intent(in) :: step,broad 4853 class(ebands_t),intent(in) :: ebands 4854 type(crystal_t),intent(in) :: cryst 4855 4856 !Local variables------------------------------- 4857 !scalars 4858 integer :: ik_ibz,ibc,ibv,spin,nw,nband_k,nbv,nproc,my_rank,cnt,mpierr,bcorr !iw, unt, 4859 real(dp) :: wtk,wmax 4860 type(htetra_t) :: tetra 4861 character(len=500) :: msg 4862 !character(len=fnlen) :: path 4863 !arrays 4864 integer :: val_idx(ebands%nkpt,ebands%nsppol) 4865 real(dp) :: eminmax(2,ebands%nsppol) 4866 real(dp),allocatable :: cvmw(:),wdt(:,:) 4867 4868 ! ********************************************************************* 4869 4870 ierr = 0 4871 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 4872 4873 jdos%nsppol = ebands%nsppol 4874 jdos%nkibz = ebands%nkpt 4875 4876 ! Find the valence band index for each k and spin. 4877 val_idx = ebands_get_valence_idx(ebands) 4878 4879 do spin=1,ebands%nsppol 4880 if (any(val_idx(:,spin) /= val_idx(1,spin))) then 4881 write(msg,'(a,i0,a)')& 4882 'Trying to compute JDOS with a metallic band structure for spin: ',spin,& 4883 'Assuming you know what you are doing, continuing anyway! ' 4884 ABI_COMMENT(msg) 4885 end if 4886 end do 4887 4888 ! Compute the linear mesh so that it encloses all bands. 4889 !if (.not. present(mesh)) then 4890 eminmax = ebands_get_minmax(ebands, "eig") 4891 wmax = maxval(eminmax(2,:) - eminmax(1,:)) 4892 nw = nint(wmax/step) + 1 4893 ABI_MALLOC(jdos%mesh, (nw)) 4894 jdos%mesh = arth(zero, step, nw) 4895 4896 jdos%nw = nw 4897 jdos%intmeth = intmeth 4898 jdos%broad = broad 4899 4900 !if (ebands%nkpt == 1) then 4901 ! ABI_COMMENT("Cannot use tetrahedra for e-DOS when nkpt == 1. Switching to gaussian method") 4902 ! jdos%intmeth = 1 4903 !end if 4904 4905 ABI_CALLOC(jdos%values, (nw, ebands%nsppol)) 4906 4907 select case (intmeth) 4908 case (1) 4909 ! Gaussian 4910 ABI_MALLOC(cvmw, (nw)) 4911 4912 cnt = 0 4913 do spin=1,ebands%nsppol 4914 do ik_ibz=1,ebands%nkpt 4915 wtk = ebands%wtk(ik_ibz) 4916 nband_k = ebands%nband(ik_ibz + (spin-1)*ebands%nkpt) 4917 nbv = val_idx(ik_ibz, spin) 4918 4919 do ibv=1,nbv 4920 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle 4921 do ibc=nbv+1,nband_k 4922 cvmw = ebands%eig(ibc,ik_ibz,spin) - ebands%eig(ibv,ik_ibz,spin) - jdos%mesh 4923 jdos%values(:, spin) = jdos%values(:, spin) + wtk * gaussian(cvmw, broad) 4924 end do 4925 end do 4926 4927 end do ! ik_ibz 4928 end do ! spin 4929 4930 ABI_FREE(cvmw) 4931 call xmpi_sum(jdos%values, comm, mpierr) 4932 4933 case (2, -2) 4934 ! Tetrahedron method 4935 if (any(ebands%nband /= ebands%nband(1)) ) then 4936 ABI_WARNING('For tetrahedra, nband(:) must be constant') 4937 ierr = ierr + 1 4938 end if 4939 if (ierr/=0) return 4940 4941 tetra = tetra_from_kptrlatt(cryst, ebands%kptopt, ebands%kptrlatt, & 4942 ebands%nshiftk, ebands%shiftk, ebands%nkpt, ebands%kptns, comm, msg, ierr) 4943 if (ierr /= 0) then 4944 call tetra%free(); return 4945 end if 4946 4947 ! For each spin and band, interpolate over kpoints, 4948 ! calculate integration weights and DOS contribution. 4949 ABI_MALLOC(cvmw, (jdos%nkibz)) 4950 ABI_MALLOC(wdt, (nw, 2)) 4951 4952 bcorr = 0; if (intmeth == -2) bcorr = 1 4953 cnt = 0 4954 do spin=1,ebands%nsppol 4955 nbv = val_idx(1, spin) 4956 do ibv=1,nbv 4957 do ibc=nbv+1,ebands%mband 4958 ! For each (c, v) get its contribution 4959 cvmw = ebands%eig(ibc,:,spin) - ebands%eig(ibv,:,spin) 4960 do ik_ibz=1,ebands%nkpt 4961 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! mpi-parallelism 4962 4963 ! Calculate integration weights at each irred k-point (Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]) 4964 call tetra%get_onewk(ik_ibz, bcorr, nw, ebands%nkpt, cvmw, jdos%mesh(0), jdos%mesh(nw), one, wdt) 4965 jdos%values(:,spin) = jdos%values(:,spin) + wdt(:, 1) * ebands%wtk(ik_ibz) 4966 end do 4967 end do ! ibc 4968 end do ! ibv 4969 end do ! spin 4970 4971 call xmpi_sum(jdos%values, comm, mpierr) 4972 4973 ! Free memory 4974 ABI_FREE(wdt) 4975 ABI_FREE(cvmw) 4976 call tetra%free() 4977 4978 case default 4979 ABI_ERROR(sjoin("Wrong integration method:", itoa(intmeth))) 4980 end select 4981 4982 if (ebands%nsppol == 1) then 4983 jdos%values(0,:) = two * jdos%values(1,:) 4984 else 4985 jdos%values(0,:) = sum(jdos%values(1:2, :), dim=2) 4986 end if 4987 4988 ! Write data. 4989 !if (my_rank == 0) then 4990 ! path = "jdos_gauss.data"; if (intmeth == 2) path = "jdos_tetra.data" 4991 ! if (open_file(path, msg, newunit=unt, form="formatted", action="write") /= 0) then 4992 ! ABI_ERROR(msg) 4993 ! end if 4994 ! do iw=1,nw 4995 ! write(unt,*)jdos%mesh(iw),(jdos(iw,spin), spin=1,ebands%nsppol) 4996 ! end do 4997 ! close(unt) 4998 !end if 4999 5000 end function ebands_get_jdos
m_ebands/ebands_get_minmax [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_minmax
FUNCTION
Report the min and max value over k-points and bands of (eig|occ|doccde) for each spin. Cannot use F90 array syntax due to the internal storage used in abinit.
INPUTS
ebands<ebands_t>=The object describing the band structure. arr_name=The name of the array whose min and Max value has to be calculated. Possible values: 'occ', 'eig' 'doccde'
OUTPUT
minmax(2,ebands%nsppol)=For each spin the min and max value of the quantity specified by "arr_name"
SOURCE
2117 function ebands_get_minmax(ebands, arr_name) result(minmax) 2118 2119 !Arguments ------------------------------------ 2120 !scalars 2121 class(ebands_t),target,intent(in) :: ebands 2122 character(len=*),intent(in) :: arr_name 2123 !arrays 2124 real(dp) :: minmax(2,ebands%nsppol) 2125 2126 !Local variables------------------------------- 2127 !scalars 2128 integer :: band,ikpt,spin,nband_k 2129 real(dp) :: datum 2130 !arrays 2131 real(dp), ABI_CONTIGUOUS pointer :: rdata(:,:,:) 2132 2133 ! ************************************************************************* 2134 2135 select case (tolower(arr_name)) 2136 case ('occ') 2137 rdata => ebands%occ 2138 case ('eig') 2139 rdata => ebands%eig 2140 case ('doccde') 2141 rdata => ebands%doccde 2142 case default 2143 ABI_BUG(sjoin('Wrong arr_name:', arr_name)) 2144 end select 2145 2146 minmax(1,:)=greatest_real 2147 minmax(2,:)=smallest_real 2148 2149 do spin=1,ebands%nsppol 2150 do ikpt=1,ebands%nkpt 2151 nband_k=ebands%nband(ikpt+(spin-1)*ebands%nkpt) 2152 do band=1,nband_k 2153 datum=rdata(band,ikpt,spin) 2154 minmax(1,spin)=MIN(minmax(1,spin),datum) 2155 minmax(2,spin)=MAX(minmax(2,spin),datum) 2156 end do 2157 end do 2158 end do 2159 2160 end function ebands_get_minmax
m_ebands/ebands_get_muT_with_fd [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_muT_with_fd
FUNCTION
Compute the Fermi level for different temperatures using Fermi-Dirac occupation function (physical T) Use ebands%nelect provided in input. Does not change input ebands.
INPUTS
OUTPUT
SOURCE
2669 subroutine ebands_get_muT_with_fd(self, ntemp, kTmesh, spinmagntarget, prtvol, mu_e, comm) 2670 2671 !Arguments ------------------------------------ 2672 !scalars 2673 class(ebands_t),intent(in) :: self 2674 integer,intent(in) :: ntemp, prtvol, comm 2675 real(dp),intent(in) :: spinmagntarget 2676 real(dp),intent(in) :: kTmesh(ntemp) 2677 real(dp),intent(out) :: mu_e(ntemp) 2678 2679 !Local variables------------------------------- 2680 !scalars 2681 integer,parameter :: occopt3 = 3 2682 integer :: ierr, it, nprocs, my_rank 2683 real(dp) :: nelect, cpu, wall, gflops 2684 type(ebands_t) :: tmp_ebands 2685 character(len=500) :: msg 2686 2687 ! ************************************************************************* 2688 2689 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 2690 call cwtime(cpu, wall, gflops, "start") 2691 2692 call ebands_copy(self, tmp_ebands) 2693 2694 mu_e = zero 2695 2696 do it=1,ntemp 2697 if (mod(it, nprocs) /= my_rank) cycle ! MPI parallelism inside comm. 2698 2699 ! Use Fermi-Dirac occopt 2700 call ebands_set_scheme(tmp_ebands, occopt3, kTmesh(it), spinmagntarget, prtvol) 2701 mu_e(it) = tmp_ebands%fermie 2702 ! 2703 ! Check that the total number of electrons is correct 2704 ! This is to trigger problems as the routines that calculate the occupations in ebands_set_extrael 2705 ! are different from the occ_fd that will be used in the rest of the code. 2706 nelect = ebands_calc_nelect(tmp_ebands, kTmesh(it), mu_e(it)) 2707 2708 if (abs(nelect - self%nelect) > tol6) then 2709 ! For T = 0 the number of occupied states goes in discrete steps (according to the k-point sampling) 2710 ! for finite doping its hard to find nelect that exactly matches self%nelect. 2711 ! in this case we print a warning 2712 write(msg,'(2(a,f10.6,a), f10.6)') & 2713 'Calculated number of electrons nelect : ',nelect, ch10, & 2714 ' does not correspond with ebands%nelect: ',tmp_ebands%nelect,' for kT: ', kTmesh(it) 2715 2716 if (kTmesh(it) == zero) then 2717 ABI_WARNING(msg) 2718 else 2719 ABI_ERROR(msg) 2720 end if 2721 end if 2722 end do ! it 2723 2724 call ebands_free(tmp_ebands) 2725 call xmpi_sum(mu_e, comm, ierr) 2726 2727 call cwtime_report(" ebands_get_muT_with_fd", cpu, wall, gflops, end_str=ch10) 2728 2729 end subroutine ebands_get_muT_with_fd
m_ebands/ebands_get_occupied [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_occupied
FUNCTION
For each k-point and spin polarisation, report the band index after which the occupation numbers are less than tol_occ.
INPUTS
ebands<ebands_t>=The object describing the band structure. tol_occ[Optional]=Tollerance on the occupation factors.
OUTPUT
NOTES
We assume that the occupation factors monotonically decrease as a function of energy. This is not always true for every smearing technique implemented in Abinit. CP: this also not true for occopt 9
SOURCE
1799 pure function ebands_get_occupied(ebands, tol_occ) result(occ_idx) 1800 1801 !Arguments ------------------------------------ 1802 !scalars 1803 real(dp),optional,intent(in) :: tol_occ 1804 class(ebands_t),intent(in) :: ebands 1805 !arrays 1806 integer :: occ_idx(ebands%nkpt,ebands%nsppol) 1807 1808 !Local variables------------------------------- 1809 integer :: band,ikpt,spin,idx,nband_k 1810 real(dp) :: tol_ 1811 1812 ! ************************************************************************* 1813 1814 tol_=tol8; if (PRESENT(tol_occ)) tol_=tol_occ 1815 1816 do spin=1,ebands%nsppol 1817 do ikpt=1,ebands%nkpt 1818 nband_k = ebands%nband(ikpt+(spin-1)*ebands%nkpt) 1819 1820 idx=0 1821 do band=1,nband_k 1822 if (ebands%occ(band,ikpt,spin) < ABS(tol_)) then 1823 idx=band; EXIT 1824 end if 1825 end do 1826 occ_idx(ikpt,spin)=idx-1 1827 if (idx==1) occ_idx(ikpt,spin)=idx 1828 if (idx==0) occ_idx(ikpt,spin)=nband_k 1829 1830 end do 1831 end do 1832 1833 end function ebands_get_occupied
m_ebands/ebands_get_valence_idx [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_get_valence_idx
FUNCTION
For each k-point and spin polarisation, report:
1) the index of the valence in case of semiconductors at T = 0
2) (band_k - 1) where band_k is the first band whose energy is > Fermi energy + told_fermi
using the value of the Fermi level.
INPUTS
ebands<ebands_t>=The object describing the band structure. tol_fermi[optional]
OUTPUT
SOURCE
1526 pure function ebands_get_valence_idx(ebands, tol_fermi) result(val_idx) 1527 1528 !Arguments ------------------------------------ 1529 !scalars 1530 real(dp),optional,intent(in) :: tol_fermi 1531 class(ebands_t),intent(in) :: ebands 1532 !arrays 1533 integer :: val_idx(ebands%nkpt,ebands%nsppol) 1534 1535 !Local variables------------------------------- 1536 integer :: band,ikpt,spin,idx,nband_k 1537 real(dp) :: tol_ 1538 1539 ! ************************************************************************* 1540 1541 tol_ = tol6; if (present(tol_fermi)) tol_ = tol_fermi 1542 1543 do spin=1,ebands%nsppol 1544 do ikpt=1,ebands%nkpt 1545 if (ebands%occopt == 9) then 1546 val_idx(ikpt,spin) = ebands%ivalence 1547 else 1548 nband_k = ebands%nband(ikpt+(spin-1)*ebands%nkpt) 1549 idx = 0 1550 do band=1,nband_k 1551 if (ebands%eig(band,ikpt,spin) > ebands%fermie + abs(tol_)) then 1552 idx = band; exit 1553 end if 1554 end do 1555 val_idx(ikpt,spin) = idx - 1 1556 if (idx == 1) val_idx(ikpt, spin) = idx 1557 if (idx == 0) val_idx(ikpt, spin) = nband_k 1558 end if 1559 end do 1560 end do 1561 1562 end function ebands_get_valence_idx
m_ebands/ebands_has_metal_scheme [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_metallic_scheme
FUNCTION
Returns .TRUE. if metallic occupation scheme is used. Note that this does not imply that the system is metallic.
INPUTS
ebands<ebands_t>=The ebands_t datatype
SOURCE
2178 pure logical function ebands_has_metal_scheme(ebands) result(ans) 2179 2180 !Arguments ------------------------------------ 2181 class(ebands_t),intent(in) :: ebands 2182 2183 ! ************************************************************************* 2184 2185 ans = (any(ebands%occopt == [3, 4, 5, 6, 7, 8, 9])) 2186 2187 end function ebands_has_metal_scheme
m_ebands/ebands_init [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_init
FUNCTION
This subroutine initializes the ebands_t structured datatype
INPUTS
bantot=total number of bands (=sum(nband(:)) doccde(bantot)=derivative of the occupation numbers with respect to the energy (Ha) eig(bantot)=eigenvalues (hartree) istwfk(nkpt)=parameter that describes the storage of wfs. ivalence = index of the valence band separating thermalized excited holes from excited thermalized excited electrons kptns(3,nkpt)=k points in terms of recip primitive translations nband(nkpt*nsppol)=number of bands nelect=Number of electrons. ne_qFD, nh_qFD= number of electrons (holes resp.) excited in the bands above band index ivalence + 1 (below ivalence resp.) nkpt=number of k points npwarr(nkpt)=number of planewaves at each k point nsppol=1 for unpolarized, 2 for spin-polarized nspinor=Number of spinor components occopt=Occupation options (see input variable) occ(bantot)=occupation numbers tphysel=Physical temperature (input variable) tsmear=Temperature of smearing. wtk(nkpt)=weight assigned to each k point cellcharge=Additional charge added to the unit cell (input variable). kptopt=Option for k-point generation (see input variable) kptrlatt_orig=Original value of kptrlatt given in input nshiftk_orig=Original number of shifts given in input shiftk_orig(3,nshiftk_orig)=Original set of shifts given in input kptrlatt=Value of kptrlatt after inkpts nshiftk=Number of shifts after inkpts shiftk(3,nshiftk)=Set of shifts after inkpts.
OUTPUT
ebands<ebands_t>=the ebands_t datatype
SOURCE
799 subroutine ebands_init(bantot, ebands, nelect, ne_qFD, nh_qFD, ivalence, doccde, eig, istwfk, kptns, & 800 nband, nkpt, npwarr, nsppol, nspinor, tphysel, tsmear, occopt, occ, wtk, & 801 cellcharge, kptopt, kptrlatt_orig, nshiftk_orig, shiftk_orig, kptrlatt, nshiftk, shiftk) 802 803 !Arguments ------------------------------------ 804 !scalars 805 integer,intent(in) :: bantot,nkpt,nsppol,nspinor,occopt,ivalence 806 real(dp),intent(in) :: nelect,ne_qFD,nh_qFD,tphysel,tsmear 807 type(ebands_t),intent(out) :: ebands 808 !arrays 809 integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol),npwarr(nkpt) 810 real(dp),intent(in) :: doccde(bantot),eig(bantot),kptns(3,nkpt),occ(bantot) 811 real(dp),intent(in) :: wtk(nkpt) 812 integer,intent(in) :: kptopt, nshiftk_orig, nshiftk 813 real(dp),intent(in) :: cellcharge 814 integer,intent(in) :: kptrlatt_orig(3,3),kptrlatt(3,3) 815 real(dp),intent(in) :: shiftk_orig(3,nshiftk_orig),shiftk(3,nshiftk) 816 817 ! ************************************************************************* 818 819 ! Copy the scalars 820 ! MG TODO here there is a inconsistency in the way occ are treated in the header 821 ! (only the states used, bantot. are saved, and the way occ. and energies 822 ! are passed to routines (mband,nkpt,nsppol). It might happen that bantot<mband*nktp*nsppol 823 ! this should not lead to problems since arrays are passed by reference 824 ! anyway the treatment of these arrays have to be rationalized 825 ebands%bantot = bantot 826 ebands%mband = MAXVAL(nband(1:nkpt*nsppol)) 827 ebands%nkpt = nkpt 828 ebands%nspinor= nspinor 829 ebands%nsppol = nsppol 830 ebands%occopt = occopt 831 832 ebands%entropy= zero 833 ebands%fermie = zero 834 ebands%fermih = zero 835 ebands%ivalence =ivalence 836 ebands%nelect = nelect 837 ebands%ne_qFD =ne_qFD 838 ebands%nh_qFD =nh_qFD 839 ebands%tphysel= tphysel 840 ebands%tsmear = tsmear 841 842 ! Allocate the components 843 ABI_MALLOC(ebands%nband, (nkpt*nsppol)) 844 ABI_MALLOC(ebands%istwfk, (nkpt)) 845 ABI_MALLOC(ebands%npwarr, (nkpt)) 846 ABI_MALLOC(ebands%kptns, (3, nkpt)) 847 848 ! Copy the arrays 849 ebands%nband(1:nkpt*nsppol) = nband(1:nkpt*nsppol) 850 ebands%istwfk(1:nkpt) = istwfk(1:nkpt) 851 ebands%npwarr(1:nkpt) = npwarr(1:nkpt) 852 ebands%kptns(1:3,1:nkpt) = kptns(1:3,1:nkpt) 853 854 ! In ebands, energies and occupations are stored in a matrix (mband,nkpt,nsppol). 855 ! put_eneocc_vect is used to reshape the values stored in vectorial form. 856 ABI_MALLOC(ebands%eig , (ebands%mband, nkpt, nsppol)) 857 ABI_MALLOC(ebands%occ , (ebands%mband, nkpt, nsppol)) 858 ABI_MALLOC(ebands%doccde, (ebands%mband, nkpt, nsppol)) 859 860 call put_eneocc_vect(ebands,'eig', eig ) 861 call put_eneocc_vect(ebands,'occ', occ ) 862 call put_eneocc_vect(ebands,'doccde',doccde) 863 864 ABI_MALLOC(ebands%wtk, (nkpt)) 865 ebands%wtk(1:nkpt) = wtk(1:nkpt) 866 867 ebands%kptopt = kptopt 868 ebands%nshiftk_orig = nshiftk_orig 869 ebands%nshiftk = nshiftk 870 ebands%cellcharge = cellcharge 871 ebands%kptrlatt_orig = kptrlatt_orig 872 ebands%kptrlatt = kptrlatt 873 874 call alloc_copy(shiftk_orig, ebands%shiftk_orig) 875 call alloc_copy(shiftk, ebands%shiftk) 876 877 end subroutine ebands_init
m_ebands/ebands_interp_kmesh [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_interp_kmesh
FUNCTION
Interpolate energies on a k-mesh.
INPUTS
ebands<ebands_t> = Object with input energies. cryst<crystal_t> = Crystalline structure. params(:): params(0): interpolation type. 1 for star-functions 2 for star-functions with group velocities if star-functions: params(2): Ratio between star functions and ab-initio k-points. params(3:4): Activate Fourier filtering (Eq 9 of PhysRevB.61.1639) if params(2) > tol6 params(3)=rcut, params(4) = rsigma intp_kptrlatt(3,3) = New k-mesh intp_nshiftk= Number of shifts in new k-mesh. intp_shiftk(3,intp_nshiftk) = Shifts in new k-mesh. band_block(2)=Initial and final band index. If [0,0], all bands are used This is a global variable i.e. all MPI procs must call the routine with the same value. comm=MPI communicator
OUTPUT
New ebands_t object with interpolated energies.
NOTES
Fermi level and occupation factors of the interpolate bands are not recomputed by this routine. This operation is delegated to the caller.
SOURCE
4208 type(ebands_t) function ebands_interp_kmesh(ebands, cryst, params, intp_kptrlatt, intp_nshiftk, intp_shiftk, & 4209 band_block, comm, out_prefix) result(new) 4210 4211 !Arguments ------------------------------------ 4212 !scalars 4213 integer,intent(in) :: intp_nshiftk,comm 4214 type(ebands_t),intent(in) :: ebands 4215 type(crystal_t),intent(in) :: cryst 4216 character(len=*),optional,intent(in) :: out_prefix 4217 !arrays 4218 integer,intent(in) :: intp_kptrlatt(3,3),band_block(2) 4219 real(dp),intent(in) :: params(:) 4220 real(dp),intent(in) :: intp_shiftk(3,intp_nshiftk) 4221 4222 !Local variables------------------------------- 4223 !scalars 4224 integer,parameter :: master = 0 4225 integer :: ik_ibz,spin,new_bantot,new_mband,cplex,itype,nb,ib 4226 integer :: nprocs,my_rank,cnt,ierr,band,new_nkbz,new_nkibz,new_nshiftk 4227 integer :: ncid 4228 type(skw_t) :: skw 4229 !arrays 4230 integer :: new_kptrlatt(3,3),my_bblock(2) 4231 integer,allocatable :: new_istwfk(:),new_nband(:,:),new_npwarr(:) 4232 real(dp),allocatable :: new_shiftk(:,:),new_kibz(:,:),new_kbz(:,:),new_wtk(:) 4233 real(dp),allocatable :: new_doccde(:),new_eig(:),new_occ(:) 4234 4235 ! ********************************************************************* 4236 4237 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 4238 itype = nint(params(1)) 4239 my_bblock = band_block; if (all(band_block == 0)) my_bblock = [1, ebands%mband] 4240 nb = my_bblock(2) - my_bblock(1) + 1 4241 4242 ! Get ibz, new shifts and new kptrlatt. 4243 call kpts_ibz_from_kptrlatt(cryst, intp_kptrlatt, ebands%kptopt, intp_nshiftk, intp_shiftk, & 4244 new_nkibz, new_kibz, new_wtk, new_nkbz, new_kbz, new_kptrlatt=new_kptrlatt, new_shiftk=new_shiftk) 4245 new_nshiftk = size(new_shiftk, dim=2) 4246 4247 ! Initialize new ebands_t in new IBZ 4248 ABI_MALLOC(new_istwfk, (new_nkibz)) 4249 new_istwfk = 1 4250 !do ik_ibz=1,new_nkibz 4251 ! new_istwfk(ik_ibz) = set_istwfk(new%kptns(:, ik_ibz)) 4252 !end do 4253 ABI_MALLOC(new_nband, (new_nkibz, ebands%nsppol)) 4254 new_nband = nb 4255 ABI_MALLOC(new_npwarr, (new_nkibz)) 4256 new_npwarr = maxval(ebands%npwarr) 4257 new_bantot = sum(new_nband); new_mband = maxval(new_nband) 4258 ABI_CALLOC(new_doccde, (new_bantot)) 4259 ABI_CALLOC(new_eig, (new_bantot)) 4260 ABI_CALLOC(new_occ, (new_bantot)) 4261 4262 call ebands_init(new_bantot, new, ebands%nelect, ebands%ne_qFD,ebands%nh_qFD,ebands%ivalence,& 4263 new_doccde, new_eig, new_istwfk, new_kibz,& 4264 new_nband, new_nkibz, new_npwarr, ebands%nsppol, ebands%nspinor, ebands%tphysel, ebands%tsmear,& 4265 ebands%occopt, new_occ, new_wtk, & 4266 ebands%cellcharge, ebands%kptopt, intp_kptrlatt, intp_nshiftk, intp_shiftk, new_kptrlatt, new_nshiftk, new_shiftk) 4267 4268 ! Get fermi level from input ebands. 4269 new%fermie = ebands%fermie 4270 new%fermih = ebands%fermih 4271 4272 ABI_FREE(new_kibz) 4273 ABI_FREE(new_wtk) 4274 ABI_FREE(new_shiftk) 4275 ABI_FREE(new_kbz) 4276 ABI_FREE(new_istwfk) 4277 ABI_FREE(new_nband) 4278 ABI_FREE(new_npwarr) 4279 ABI_FREE(new_doccde) 4280 ABI_FREE(new_eig) 4281 ABI_FREE(new_occ) 4282 4283 ! Build SKW object for all bands. 4284 if (itype == 1 .or. itype == 2) then 4285 cplex = 1; if (kpts_timrev_from_kptopt(ebands%kptopt) == 0) cplex = 2 4286 skw = skw_new(cryst, params(2:), cplex, ebands%mband, ebands%nkpt, ebands%nsppol, ebands%kptns, ebands%eig, my_bblock, comm) 4287 !if (itype == 2) then 4288 ! ABI_CALLOC(new%velocity,(3,new%mband,new%nkpt,new%nsppol)) 4289 !end if 4290 else 4291 ABI_ERROR(sjoin("Wrong einterp params(1):", itoa(itype))) 4292 end if 4293 4294 ! Interpolate eigenvalues and velocities. 4295 new%eig = zero; cnt = 0 4296 do spin=1,new%nsppol 4297 do ik_ibz=1,new%nkpt 4298 do ib=1,nb 4299 cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism. 4300 ! Note the difference between band and ib index if band_block. 4301 band = my_bblock(1) + ib - 1 4302 select case (itype) 4303 case (1) 4304 call skw%eval_bks(band, new%kptns(:,ik_ibz), spin, new%eig(ib,ik_ibz,spin)) 4305 !case (2) 4306 ! call skw%eval_bks(band, new%kptns(:,ik_ibz), spin, new%eig(ib,ik_ibz,spin), new%velocity(:,ib,ik_ibz,spin)) 4307 case default 4308 ABI_ERROR(sjoin("Wrong params(1):", itoa(itype))) 4309 end select 4310 end do 4311 end do 4312 end do 4313 call xmpi_sum(new%eig, comm, ierr) 4314 !if (itype == 2) call xmpi_sum(new%velocity, comm, ierr) 4315 4316 ! Sort eigvalues_k in ascending order to be compatible with other ebands routines. 4317 call ebands_sort(new) 4318 !call ebands_update_occ(new, dtset%spinmagntarget, prtvol=dtset%prtvol) 4319 4320 if (my_rank == master .and. itype == 1 .and. present(out_prefix)) then 4321 ! Write ESKW file with crystal and (interpolated) band structure energies. 4322 !call wrtout(ab_out, sjoin("- Writing interpolated bands to file:", strcat(prefix, tag))) 4323 ! Write crystal and (interpolated) band structure energies. 4324 NCF_CHECK(nctk_open_create(ncid, strcat(out_prefix, "_ESKW.nc"), xmpi_comm_self)) 4325 NCF_CHECK(cryst%ncwrite(ncid)) 4326 NCF_CHECK(ebands_ncwrite(new, ncid)) 4327 4328 ! TODO 4329 !NCF_CHECK(skw%ncwrite(ncid)) 4330 ! Define variables specific to SKW algo. 4331 !ncerr = nctk_def_arrays(ncid, [ & 4332 ! nctkarr_t("band_block", "int", "two"), & 4333 ! nctkarr_t("einterp", "dp", "four")], defmode=.True.) 4334 !NCF_CHECK(ncerr) 4335 4336 ! Write data. 4337 !NCF_CHECK(nctk_set_datamode(ncid)) 4338 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "band_block"), band_block)) 4339 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "einterp"), params)) 4340 NCF_CHECK(nf90_close(ncid)) 4341 end if 4342 4343 call skw%free() 4344 4345 end function ebands_interp_kmesh
m_ebands/ebands_interp_kpath [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_interp_kpath
FUNCTION
Interpolate energies on a k-path
INPUTS
ebands<ebands_t> = Object with input energies. cryst<crystal_t> = Crystalline structure. kpath<kpath_t> = Object describing the k-path params(:): params(1): 1 for SKW. band_block(2)=Initial and final band index to be interpolated. [0,0] if all bands are used. This is a global variable i.e. all MPI procs must call the routine with the same value. comm=MPI communicator
OUTPUT
New ebands_t object with interpolated energies.
SOURCE
4372 type(ebands_t) function ebands_interp_kpath(ebands, cryst, kpath, params, band_block, comm) result(new) 4373 4374 !Arguments ------------------------------------ 4375 !scalars 4376 integer,intent(in) :: comm 4377 type(ebands_t),intent(in) :: ebands 4378 type(crystal_t),intent(in) :: cryst 4379 type(kpath_t),intent(in) :: kpath 4380 !arrays 4381 integer,intent(in) :: band_block(2) 4382 real(dp),intent(in) :: params(:) 4383 4384 !Local variables------------------------------- 4385 !scalars 4386 integer,parameter :: new_nshiftk=1 4387 integer :: ik_ibz,spin,new_bantot,new_mband,cplex 4388 integer :: nprocs,my_rank,cnt,ierr,band,new_nkibz,itype,nb,ib, new_kptopt 4389 type(skw_t) :: skw 4390 !arrays 4391 integer,parameter :: new_kptrlatt(3,3)=0 4392 integer :: my_bblock(2) 4393 integer,allocatable :: new_istwfk(:),new_nband(:,:),new_npwarr(:) 4394 real(dp),parameter :: new_shiftk(3,1) = zero 4395 real(dp),allocatable :: new_wtk(:),new_doccde(:),new_eig(:),new_occ(:) 4396 4397 ! ********************************************************************* 4398 4399 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 4400 itype = nint(params(1)) 4401 my_bblock = band_block; if (all(band_block == 0)) my_bblock = [1, ebands%mband] 4402 nb = my_bblock(2) - my_bblock(1) + 1 4403 4404 if (ebands%nkpt == 1) then 4405 ABI_WARNING("Cannot interpolate band energies when nkpt = 1. Returning") 4406 return 4407 end if 4408 4409 ! Initialize new ebands_t. 4410 new_nkibz = kpath%npts 4411 ABI_MALLOC(new_istwfk, (new_nkibz)) 4412 new_istwfk = 1 4413 ABI_MALLOC(new_nband, (new_nkibz, ebands%nsppol)) 4414 new_nband = nb 4415 ABI_MALLOC(new_npwarr, (new_nkibz)) 4416 new_npwarr = maxval(ebands%npwarr) 4417 new_bantot = sum(new_nband); new_mband = maxval(new_nband) 4418 ABI_CALLOC(new_eig, (new_bantot)) 4419 ABI_CALLOC(new_doccde, (new_bantot)) 4420 ABI_CALLOC(new_occ, (new_bantot)) 4421 ABI_CALLOC(new_wtk, (new_nkibz)) 4422 4423 ! Needed by AbiPy to understand that we have a k-path instead of a mesh. 4424 new_kptopt = -kpath%nbounds 4425 4426 call ebands_init(new_bantot, new, ebands%nelect, ebands%ne_qFD,ebands%nh_qFD,ebands%ivalence, & 4427 new_doccde, new_eig, new_istwfk, kpath%points, & 4428 new_nband, new_nkibz, new_npwarr, ebands%nsppol, ebands%nspinor, ebands%tphysel, ebands%tsmear, & 4429 ebands%occopt, new_occ, new_wtk,& 4430 ebands%cellcharge, new_kptopt, new_kptrlatt, new_nshiftk, new_shiftk, new_kptrlatt, new_nshiftk, new_shiftk) 4431 4432 new%fermie = ebands%fermie 4433 new%fermih = ebands%fermih 4434 4435 ABI_FREE(new_wtk) 4436 ABI_FREE(new_istwfk) 4437 ABI_FREE(new_nband) 4438 ABI_FREE(new_npwarr) 4439 ABI_FREE(new_doccde) 4440 ABI_FREE(new_eig) 4441 ABI_FREE(new_occ) 4442 4443 ! Build SKW object for all bands. 4444 select case (itype) 4445 case (1) 4446 cplex = 1; if (kpts_timrev_from_kptopt(ebands%kptopt) == 0) cplex = 2 4447 skw = skw_new(cryst, params(2:), cplex, ebands%mband, ebands%nkpt, ebands%nsppol, ebands%kptns, ebands%eig, & 4448 my_bblock, comm) 4449 4450 case default 4451 ABI_ERROR(sjoin("Wrong einterp params(1):", itoa(itype))) 4452 end select 4453 4454 ! Interpolate eigenvalues. 4455 new%eig = zero; cnt = 0 4456 do spin=1,new%nsppol 4457 do ik_ibz=1,new%nkpt 4458 do ib=1,nb 4459 cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism. 4460 ! Note the difference between band and ib index if band_block. 4461 band = my_bblock(1) + ib - 1 4462 select case (itype) 4463 case (1) 4464 call skw%eval_bks(band, new%kptns(:,ik_ibz), spin, new%eig(ib,ik_ibz,spin)) 4465 case default 4466 ABI_ERROR(sjoin("Wrong einterp params(1):", itoa(itype))) 4467 end select 4468 end do 4469 end do 4470 end do 4471 call xmpi_sum(new%eig, comm, ierr) 4472 4473 ! Sort eigvalues_k in ascending order to be compatible with other ebands routines. 4474 call ebands_sort(new) 4475 4476 call skw%free() 4477 4478 end function ebands_interp_kpath
m_ebands/ebands_interpolate_kpath [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_interpolate_kpath
FUNCTION
INPUTS
dtset<dataset_type>=Abinit dataset band_block(2)=Initial and final band index to be interpolated. [0,0] if all bands are used. This is a global variable i.e. all MPI procs must call the routine with the same value.
OUTPUT
SOURCE
5824 subroutine ebands_interpolate_kpath(ebands, dtset, cryst, band_block, prefix, comm) 5825 5826 !Arguments ------------------------------------ 5827 !scalars 5828 type(ebands_t),intent(in) :: ebands 5829 type(dataset_type),intent(in) :: dtset 5830 type(crystal_t),intent(in) :: cryst 5831 integer,intent(in) :: comm 5832 character(len=*),intent(in) :: prefix 5833 !arrays 5834 integer,intent(in) :: band_block(2) 5835 5836 !Local variables------------------------------- 5837 !scalars 5838 integer,parameter :: master = 0, intp_nshiftk1 = 1 5839 integer :: my_rank, ndivsm, nbounds, itype !, spin, ik, ib, ii, jj, ierr 5840 type(ebands_t) :: ebands_kpath 5841 type(kpath_t) :: kpath 5842 character(len=500) :: tag !msg 5843 !arrays 5844 real(dp),allocatable :: bounds(:,:) 5845 5846 ! ********************************************************************* 5847 5848 my_rank = xmpi_comm_rank(comm) 5849 5850 itype = nint(dtset%einterp(1)); tag = "_SKW" 5851 tag = "_INTERP" 5852 5853 ! Generate k-path 5854 ndivsm = dtset%ndivsm 5855 if (ndivsm <= 0) then 5856 ABI_COMMENT("Setting ndivsm to 20 because variable is not given in input file") 5857 ndivsm = 20 5858 end if 5859 nbounds = dtset%nkpath 5860 if (nbounds <= 0) then 5861 ABI_COMMENT("Using hard-coded k-path because nkpath not present in input file.") 5862 nbounds = 5 5863 ABI_MALLOC(bounds, (3, 5)) 5864 bounds = reshape([zero, zero, zero, half, zero, zero, zero, half, zero, zero, zero, zero, zero, zero, half], [3,5]) 5865 else 5866 call alloc_copy(dtset%kptbounds, bounds) 5867 end if 5868 5869 kpath = kpath_new(bounds, cryst%gprimd, ndivsm) 5870 call kpath%print(header="Interpolating energies on k-path", unit=std_out) 5871 ABI_FREE(bounds) 5872 5873 ! Interpolate bands on k-path. 5874 ebands_kpath = ebands_interp_kpath(ebands, cryst, kpath, dtset%einterp, band_block, comm) 5875 5876 if (my_rank == master) then 5877 call wrtout(ab_out, sjoin("- Writing interpolated bands to file:", strcat(prefix, tag))) 5878 call ebands_write(ebands_kpath, dtset%prtebands, strcat(prefix, tag), kptbounds=kpath%bounds) 5879 end if 5880 5881 call ebands_free(ebands_kpath) 5882 call kpath%free() 5883 5884 end subroutine ebands_interpolate_kpath
m_ebands/ebands_move_alloc [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_move_alloc
FUNCTION
Transfer allocate from `from_ebands` to `to_ebands`. `from_ebands` is destroyed when the routine returns.
SOURCE
1131 subroutine ebands_move_alloc(from_ebands, to_ebands) 1132 1133 !Arguments ------------------------------------ 1134 class(ebands_t),intent(inout) :: from_ebands 1135 class(ebands_t),intent(inout) :: to_ebands 1136 1137 ! ********************************************************************* 1138 1139 call ebands_free(to_ebands) 1140 call ebands_copy(from_ebands, to_ebands) 1141 call ebands_free(from_ebands) 1142 1143 end subroutine ebands_move_alloc
m_ebands/ebands_ncwrite [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_ncwrite
FUNCTION
Writes the content of an ebands_t object to a NETCDF file according to the ETSF-IO specifications. Return nf90_noerr if success.
INPUTS
ncid =NC file handle
SOURCE
2908 integer function ebands_ncwrite(ebands, ncid) result(ncerr) 2909 2910 !Arguments ------------------------------------ 2911 !scalars 2912 integer,intent(in) :: ncid 2913 class(ebands_t),intent(in) :: ebands 2914 2915 !Local variables------------------------------- 2916 !scalars 2917 integer :: ii,nelect_int 2918 logical :: write_ngkpt 2919 character(len=etsfio_charlen) :: smearing,k_dependent 2920 !arrays 2921 integer :: ngkpt(3) 2922 2923 ! ************************************************************************* 2924 2925 smearing = nctk_string_from_occopt(ebands%occopt) 2926 2927 ! ============================================== 2928 ! === Write the dimensions specified by ETSF === 2929 ! ============================================== 2930 ncerr = nctk_def_dims(ncid, [ & 2931 nctkdim_t("max_number_of_states", ebands%mband), & 2932 nctkdim_t("number_of_spinor_components", ebands%nspinor), & 2933 nctkdim_t("number_of_spins", ebands%nsppol), & 2934 nctkdim_t("number_of_kpoints", ebands%nkpt), & 2935 nctkdim_t("nshiftk_orig", ebands%nshiftk_orig), & 2936 nctkdim_t("nshiftk", ebands%nshiftk)], & 2937 defmode=.True.) 2938 NCF_CHECK(ncerr) 2939 2940 ! FIXME 2941 ! Unofficial variables. Notes: 2942 ! 1) ETSF-IO does not support nshifts > 1 2943 ! 2) shiftk_orig, nshiftk_orig refers to the values specified in the input (most useful ones). 2944 ! 3) shiftk, kptrlatt refers to the values computed in inkpts. 2945 ! 4) Should define a protocol so that abipy understands if we have a path or a mesh. 2946 !write_kptrlatt = (SUM(ABS(ebands%kptrlatt))/=0) 2947 !write_kptrlatt = (ebands%kptopt /= 0) 2948 2949 ngkpt = 0; write_ngkpt = .False. 2950 if (isdiagmat(ebands%kptrlatt) .and. ebands%nshiftk == 1) then 2951 write_ngkpt = .True. 2952 do ii=1,3 2953 ngkpt(ii) = ebands%kptrlatt(ii, ii) 2954 end do 2955 ncerr = nctk_def_dims(ncid, nctkdim_t('ngkpt_nshiftk', ebands%nshiftk_orig)) 2956 NCF_CHECK(ncerr) 2957 end if 2958 2959 ! Define k-points 2960 ncerr = nctk_def_arrays(ncid, [& 2961 nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), & 2962 nctkarr_t("kpoint_weights", "dp", "number_of_kpoints"), & 2963 nctkarr_t("monkhorst_pack_folding", "int", "number_of_vectors") & 2964 ]) 2965 NCF_CHECK(ncerr) 2966 2967 ! Define states section. 2968 ncerr = nctk_def_arrays(ncid, [& 2969 nctkarr_t("number_of_states", "int", "number_of_kpoints, number_of_spins"), & 2970 nctkarr_t("eigenvalues", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), & 2971 nctkarr_t("occupations", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), & 2972 nctkarr_t("smearing_scheme", "char", "character_string_length") & 2973 ]) 2974 NCF_CHECK(ncerr) 2975 2976 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "number_of_electrons"]) 2977 NCF_CHECK(ncerr) 2978 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"]) 2979 NCF_CHECK(ncerr) 2980 2981 ! Some variables require the specifications of units. 2982 NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues")) 2983 NCF_CHECK(nctk_set_atomic_units(ncid, "fermi_energy")) 2984 2985 k_dependent = "no"; if (any(ebands%nband(1) /= ebands%nband)) k_dependent = "yes" 2986 NCF_CHECK(nf90_put_att(ncid, vid("number_of_states"), "k_dependent", k_dependent)) 2987 2988 ! Write data. 2989 ! 1) Electrons. 2990 ! NB: In etsf_io the number of electrons is declared as integer. 2991 ! We use abinit nelect to store the value as real(dp). 2992 nelect_int = nint(ebands%nelect) 2993 2994 NCF_CHECK(nctk_set_datamode(ncid)) 2995 NCF_CHECK(nf90_put_var(ncid, vid("fermi_energy"), ebands%fermie)) 2996 NCF_CHECK(nf90_put_var(ncid, vid("number_of_electrons"), nelect_int)) 2997 NCF_CHECK(nf90_put_var(ncid, vid("smearing_width"), ebands%tsmear)) 2998 NCF_CHECK(nf90_put_var(ncid, vid("number_of_states"), ebands%nband, count=[ebands%nkpt, ebands%nsppol])) 2999 NCF_CHECK(nf90_put_var(ncid, vid("eigenvalues"), ebands%eig)) 3000 NCF_CHECK(nf90_put_var(ncid, vid("occupations"), ebands%occ)) 3001 NCF_CHECK(nf90_put_var(ncid, vid("smearing_scheme"), smearing)) 3002 3003 ! K-points 3004 NCF_CHECK(nf90_put_var(ncid, vid("reduced_coordinates_of_kpoints"), ebands%kptns)) 3005 NCF_CHECK(nf90_put_var(ncid, vid("kpoint_weights"), ebands%wtk)) 3006 3007 if (write_ngkpt) then 3008 NCF_CHECK(nf90_put_var(ncid, vid("monkhorst_pack_folding"), ngkpt)) 3009 end if 3010 3011 ! =========================================================== 3012 ! === Write abinit-related stuff (not covered by ETSF-IO) === 3013 ! =========================================================== 3014 ! Define variables. 3015 NCF_CHECK(nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "occopt", "kptopt"], defmode=.True.)) 3016 NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "tphysel", "charge", "nelect", "extrael"])) 3017 3018 ncerr = nctk_def_arrays(ncid, nctkarr_t('istwfk', "i", 'number_of_kpoints')) 3019 NCF_CHECK(ncerr) 3020 3021 ! Abinit variables defining the K-point sampling. 3022 ncerr = nctk_def_arrays(ncid, [& 3023 nctkarr_t('kptrlatt_orig', "i", 'number_of_reduced_dimensions, number_of_reduced_dimensions'), & 3024 nctkarr_t('shiftk_orig', "dp", 'number_of_reduced_dimensions, nshiftk_orig'), & 3025 nctkarr_t('kptrlatt', "i", 'number_of_reduced_dimensions, number_of_reduced_dimensions'), & 3026 nctkarr_t('shiftk', "dp", 'number_of_reduced_dimensions, nshiftk') & 3027 ]) 3028 NCF_CHECK(ncerr) 3029 3030 if (write_ngkpt) then 3031 ncerr = nctk_def_arrays(ncid, nctkarr_t('ngkpt_shiftk', "dp", "number_of_reduced_dimensions, ngkpt_nshiftk")) 3032 NCF_CHECK(ncerr) 3033 end if 3034 3035 ! Write Abinit variables 3036 NCF_CHECK(nctk_set_datamode(ncid)) 3037 NCF_CHECK(nf90_put_var(ncid, vid("tphysel"), ebands%tphysel)) 3038 NCF_CHECK(nf90_put_var(ncid, vid("occopt"), ebands%occopt)) 3039 NCF_CHECK(nf90_put_var(ncid, vid("istwfk"), ebands%istwfk)) 3040 NCF_CHECK(nf90_put_var(ncid, vid("kptopt"), ebands%kptopt)) 3041 NCF_CHECK(nf90_put_var(ncid, vid("charge"), ebands%cellcharge)) 3042 NCF_CHECK(nf90_put_var(ncid, vid("extrael"), ebands%extrael)) 3043 NCF_CHECK(nf90_put_var(ncid, vid("nelect"), ebands%nelect)) 3044 NCF_CHECK(nf90_put_var(ncid, vid('kptrlatt_orig'), ebands%kptrlatt_orig)) 3045 NCF_CHECK(nf90_put_var(ncid, vid('shiftk_orig'), ebands%shiftk_orig)) 3046 NCF_CHECK(nf90_put_var(ncid, vid('kptrlatt'),ebands%kptrlatt)) 3047 NCF_CHECK(nf90_put_var(ncid, vid('shiftk'), ebands%shiftk)) 3048 3049 if (write_ngkpt) then 3050 NCF_CHECK(nf90_put_var(ncid, vid('ngkpt_shiftk'), ebands%shiftk_orig)) 3051 end if 3052 3053 !In the case occopt 9 3054 if (ebands%occopt == 9) then 3055 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: & 3056 "number_of_conduction_electrons", "number_of_valence_holes"]) 3057 NCF_CHECK(ncerr) 3058 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "holes_fermi_energy"]) 3059 NCF_CHECK(ncerr) 3060 NCF_CHECK(nctk_set_atomic_units(ncid,"holes_fermi_energy")) 3061 NCF_CHECK(nf90_put_var(ncid, vid("holes_fermi_energy"), ebands%fermih)) 3062 NCF_CHECK(nf90_put_var(ncid, vid("number_of_conduction_electrons"), ebands%ne_qFD)) 3063 NCF_CHECK(nf90_put_var(ncid, vid("number_of_valence_holes"), ebands%nh_qFD)) 3064 endif 3065 3066 contains 3067 integer function vid(vname) 3068 character(len=*),intent(in) :: vname 3069 vid = nctk_idname(ncid, vname) 3070 end function vid 3071 3072 end function ebands_ncwrite
m_ebands/ebands_ncwrite_path [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_ncwrite_path
FUNCTION
Writes the content of an ebands_t object to a NETCDF file
INPUTS
cryst=Crystal structure path=File name
OUTPUT
SOURCE
3092 integer function ebands_ncwrite_path(ebands, cryst, path) result(ncerr) 3093 3094 !Arguments ------------------------------------ 3095 !scalars 3096 class(ebands_t),intent(in) :: ebands 3097 type(crystal_t),intent(in) :: cryst 3098 character(len=*),intent(in) :: path 3099 3100 !Local variables------------------------------- 3101 !scalars 3102 integer :: ncid 3103 3104 ! ************************************************************************* 3105 3106 ncerr = nf90_noerr 3107 if (file_exists(path)) then 3108 NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self)) 3109 else 3110 ncerr = nctk_open_create(ncid, path, xmpi_comm_self) 3111 NCF_CHECK_MSG(ncerr, sjoin("Creating", path)) 3112 end if 3113 3114 NCF_CHECK(cryst%ncwrite(ncid)) 3115 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 3116 NCF_CHECK(nf90_close(ncid)) 3117 3118 end function ebands_ncwrite_path
m_ebands/ebands_nelect_per_spin [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_nelect_per_spin
FUNCTION
Return number of electrons in each spin channel (computed from occoputation factors if nsppol=2)
INPUTS
ebands<ebands_t>=The object describing the band structure.
OUTPUT
nelect_per_spin(ebands%nsppol)=For each spin the number of electrons (eventually fractional)
SOURCE
2068 pure function ebands_nelect_per_spin(ebands) result(nelect_per_spin) 2069 2070 !Arguments ------------------------------------ 2071 !scalars 2072 class(ebands_t),intent(in) :: ebands 2073 !arrays 2074 real(dp) :: nelect_per_spin(ebands%nsppol) 2075 2076 !Local variables------------------------------- 2077 !scalars 2078 integer :: iband,ikpt,spin 2079 2080 ! ************************************************************************* 2081 2082 nelect_per_spin = ebands%nelect 2083 if (ebands%nsppol > 1) then 2084 nelect_per_spin = zero 2085 do spin=1,ebands%nsppol 2086 do ikpt=1,ebands%nkpt 2087 do iband=1,ebands%nband(ikpt+ebands%nkpt*(spin-1)) 2088 nelect_per_spin(spin) = nelect_per_spin(spin) + ebands%wtk(ikpt) * ebands%occ(iband, ikpt, spin) 2089 end do 2090 end do 2091 end do 2092 end if 2093 2094 end function ebands_nelect_per_spin
m_ebands/ebands_print [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_print
FUNCTION
Print the content of the object.
INPUTS
ebands<ebands_t>The type containing the data. [unit]=Unit number (default: std_out) [header]=title for info [prtvol]=Verbosity level (default: 0)
OUTPUT
Only writing
SOURCE
1166 subroutine ebands_print(ebands, header, unit, prtvol) 1167 1168 !Arguments ------------------------------------ 1169 class(ebands_t),intent(in) :: ebands 1170 integer,optional,intent(in) :: prtvol, unit 1171 character(len=*),optional,intent(in) :: header 1172 1173 !Local variables------------------------------- 1174 integer :: spin, ikpt, unt, my_prtvol, ii 1175 character(len=500) :: msg 1176 ! ************************************************************************* 1177 1178 unt = std_out; if (present(unit)) unt =unit 1179 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 1180 1181 msg = ' ==== Info on the ebands_t ==== ' 1182 if (present(header)) msg=' ==== '//trim(adjustl(header))//' ==== ' 1183 call wrtout(unt, msg) 1184 1185 write(msg,'(6(a,i0,a))')& 1186 ' Number of spinorial components ...... ',ebands%nspinor,ch10,& 1187 ' Number of spin polarizations ........ ',ebands%nsppol,ch10,& 1188 ' Number of k-points in the IBZ ....... ',ebands%nkpt,ch10,& 1189 ' kptopt .............................. ',ebands%kptopt,ch10,& 1190 ' Maximum number of bands ............. ',ebands%mband,ch10,& 1191 ' Occupation option ................... ',ebands%occopt,ch10 1192 call wrtout(unt, msg) 1193 1194 write(msg,"(2a)")" kptrlatt .............. ",trim(ltoa(reshape(ebands%kptrlatt, [9]))) 1195 call wrtout(unt, msg) 1196 write(msg,"(2a)")" shiftk ................ ",trim(ltoa(reshape(ebands%shiftk, [3 * ebands%nshiftk]))) 1197 call wrtout(unt, msg) 1198 1199 write(msg,'(3(a,f14.2,a),4(a,f14.6,a))')& 1200 ' Number of valence electrons ......... ',ebands%nelect,ch10,& 1201 ' Extra cell charge (from GS run)...... ',ebands%cellcharge,ch10,& 1202 ' Extra electrons (after GS run)....... ',ebands%extrael,ch10,& 1203 ' Fermi level ........................ ',ebands%fermie,ch10,& 1204 ' Entropy ............................. ',ebands%entropy,ch10,& 1205 ' Tsmear value ........................ ',ebands%tsmear,ch10,& 1206 ' Tphysel value ....................... ',ebands%tphysel,ch10 1207 call wrtout(unt, msg) 1208 1209 if (my_prtvol > 10) then 1210 if (ebands%nsppol == 1)then 1211 call wrtout(unt, sjoin(' New occ. numbers for occopt= ', itoa(ebands%occopt),' , spin-unpolarized case.')) 1212 end if 1213 1214 do spin=1,ebands%nsppol 1215 if (ebands%nsppol == 2) then 1216 write(msg,'(a,i9,a,i0)')' New occ. numbers for occopt= ',ebands%occopt,', spin ',spin 1217 call wrtout(unt, msg) 1218 end if 1219 1220 do ikpt=1,ebands%nkpt 1221 write(msg,'(2a,i4,3a,f6.3,2a)')ch10,& 1222 ' k-point number ',ikpt,') ',trim(ktoa(ebands%kptns(:,ikpt))),'; weight: ',ebands%wtk(ikpt), ch10, & 1223 " eig (Ha), eig (eV), occ, doccde" 1224 call wrtout(unt, msg) 1225 do ii=1,ebands%nband(ikpt+(spin-1)*ebands%nkpt) 1226 write(msg,'(4(f7.3,1x))')ebands%eig(ii,ikpt,spin), ebands%eig(ii,ikpt,spin) * Ha_eV, & 1227 ebands%occ(ii,ikpt,spin), ebands%doccde(ii,ikpt,spin) 1228 call wrtout(unt, msg) 1229 end do 1230 end do !ikpt 1231 1232 end do !spin 1233 1234 end if !my_prtvol 1235 1236 end subroutine ebands_print
m_ebands/ebands_print_gaps [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_print_gaps
FUNCTION
Helper function to print gaps directly from ebands.
INPUTS
ebands<ebands_t>=Info on the band structure, the smearing technique and the physical temperature used.
OUTPUT
SOURCE
386 subroutine ebands_print_gaps(ebands, unit, header) 387 388 !Arguments ------------------------------------ 389 class(ebands_t),intent(in) :: ebands 390 integer,intent(in) :: unit 391 character(len=*),optional,intent(in) :: header 392 393 !Local variables------------------------------- 394 integer :: ierr, spin 395 type(gaps_t) :: gaps 396 397 ! ********************************************************************* 398 399 if (unit == dev_null) return 400 401 gaps = ebands_get_gaps(ebands, ierr) 402 if (ierr /= 0) then 403 do spin=1, ebands%nsppol 404 write(unit, "(a)")trim(gaps%errmsg_spin(spin)) 405 end do 406 end if 407 408 if (present(header)) then 409 call gaps%print(unit=unit, header=header) 410 else 411 call gaps%print(unit=unit) 412 end if 413 call gaps%free() 414 415 end subroutine ebands_print_gaps
m_ebands/ebands_prtbltztrp [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_prtbltztrp
FUNCTION
Output files for BoltzTraP code, which integrates Boltzmann transport quantities over the Fermi surface for different T and chemical potentials. Abinit provides all necessary input files: struct, energy, input file, and def file for the unit definitions of fortran files in BT. See http://www.icams.de/content/departments/ams/madsen/boltztrap.html
INPUTS
ebands<ebands_t>=Band structure object. cryst<cryst_t>=Info on the crystalline structure. fname_radix = radix of file names for output
OUTPUT
(only writing, printing)
SOURCE
5114 subroutine ebands_prtbltztrp(ebands, crystal, fname_radix, tau_k) 5115 5116 !Arguments ------------------------------------ 5117 !scalars 5118 type(ebands_t),intent(in) :: ebands 5119 type(crystal_t),intent(in) :: crystal 5120 character(len=fnlen), intent(in) :: fname_radix 5121 !arrays 5122 real(dp), intent(in), optional :: tau_k(ebands%nsppol,ebands%nkpt,ebands%mband) 5123 5124 !Local variables------------------------------- 5125 !scalars 5126 integer :: iout, isym, iband, isppol, ikpt, nsppol, nband 5127 real(dp),parameter :: ha2ryd=two 5128 real(dp) :: ewindow 5129 character(len=fnlen) :: filename 5130 character(len=2) :: so_suffix 5131 character(len=500) :: msg 5132 !arrays 5133 real(dp) :: nelec(ebands%nsppol) 5134 character(len=3) :: spinsuffix(ebands%nsppol) 5135 5136 ! ************************************************************************* 5137 5138 !MG FIXME The number of electrons is wrong if the file is produced in a NSCF run. 5139 ! See https://forum.abinit.org/viewtopic.php?f=19&t=3339 5140 5141 ! CP test to prevent use in case occopt = 9 5142 if (ebands%occopt==9) then 5143 write(msg,'(a)') "Boltztrap outputting not possible with occopt = 9 at the moment" 5144 ABI_ERROR(msg) 5145 end if 5146 5147 nelec = ebands_nelect_per_spin(ebands) 5148 nsppol = ebands%nsppol 5149 nband = ebands%nband(1) 5150 5151 so_suffix="" 5152 if (nsppol > 1 .or. ebands%nspinor > 1) so_suffix="so" 5153 5154 if (nsppol == 1) then 5155 spinsuffix(1) = "ns_" 5156 else 5157 spinsuffix = ["up_", "dn_"] 5158 end if 5159 5160 do isppol=1,nsppol 5161 5162 !input file for boltztrap: general info, Ef, Nelec, etc... 5163 filename= trim(fname_radix)//"_"//trim(spinsuffix(isppol))//"BLZTRP.intrans" 5164 if (open_file(filename, msg, newunit=iout, form='formatted') /= 0) then 5165 ABI_ERROR(msg) 5166 end if 5167 5168 ewindow = 1.1_dp * (ebands%fermie-minval(ebands%eig(1, :, isppol))) 5169 write (iout, '(a)') "GENE # Format of input: generic format, with Symmetries" 5170 write (iout, '(a)') "0 0 0 0.0 # iskip (not presently used) idebug setgap shiftgap" 5171 write (iout, '(E15.5,a,2F10.4,a)') ebands%fermie*ha2ryd, " 0.0005 ", ewindow*ha2ryd, nelec(isppol), & 5172 & " # Fermilevel (Ry), energy grid spacing, energy span around Fermilevel, number of electrons for this spin" 5173 write (iout, '(a)') "CALC # CALC (calculate expansion coeff), NOCALC read from file" 5174 write (iout, '(a)') "3 # lpfac, number of latt-points per k-point" 5175 write (iout, '(a)') "BOLTZ # run mode (only BOLTZ is supported)" 5176 write (iout, '(a)') ".15 # (efcut) energy range of chemical potential" 5177 write (iout, '(a)') "300. 10. # Tmax, temperature grid spacing" 5178 write (iout, '(2a)') "-1 # energyrange of bands given ",& 5179 & "individual DOS output sig_xxx and dos_xxx (xxx is band number)" 5180 write (iout, '(a)') "HISTO # DOS calculation method. Other possibility is TETRA" 5181 write (iout, '(a)') "No # not using model for relaxation time" 5182 write (iout, '(a)') "3 # Number of doping levels coefficients will be output for" 5183 write (iout, '(a)') "-1.e16 0.0d0 1.e16 # Values of doping levels (in carriers / cm^3" 5184 close(iout) 5185 5186 !files file, with association of all units for Boltztrap 5187 filename= trim(fname_radix)//"_"//trim(spinsuffix(isppol))//"BLZTRP.def" 5188 if (open_file(filename, msg, newunit=iout, form='formatted') /= 0) then 5189 ABI_ERROR(msg) 5190 end if 5191 5192 write (iout, '(3a)') "5, '", trim(fname_radix)//"_"//trim(spinsuffix(isppol))//"BLZTRP.intrans', 'old', 'formatted',0" 5193 write (iout, '(3a)') "6, '", trim(fname_radix)//"_BLZTRP", ".outputtrans', 'unknown', 'formatted',0" 5194 write (iout, '(3a)') "20,'", trim(fname_radix)//"_BLZTRP", ".struct', 'old', 'formatted',0" 5195 write (iout, '(3a)') "10,'", trim(fname_radix)//"_BLZTRP."//trim(spinsuffix(isppol))//"energy"//trim(so_suffix),& 5196 & "', 'old', 'formatted',0" 5197 if (present (tau_k)) then 5198 write (iout, '(3a)') "11,'", trim(fname_radix)//"_BLZTRP", ".tau_k', 'old', 'formatted',0" 5199 end if 5200 write (iout, '(3a)') "48,'", trim(fname_radix)//"_BLZTRP", ".engre', 'unknown', 'unformatted',0" 5201 write (iout, '(3a)') "49,'", trim(fname_radix)//"_BLZTRP", ".transdos', 'unknown', 'formatted',0" 5202 write (iout, '(3a)') "50,'", trim(fname_radix)//"_BLZTRP", ".sigxx', 'unknown', 'formatted',0" 5203 write (iout, '(3a)') "51,'", trim(fname_radix)//"_BLZTRP", ".sigxxx', 'unknown', 'formatted',0" 5204 write (iout, '(3a)') "21,'", trim(fname_radix)//"_BLZTRP", ".trace', 'unknown', 'formatted',0" 5205 write (iout, '(3a)') "22,'", trim(fname_radix)//"_BLZTRP", ".condtens', 'unknown', 'formatted',0" 5206 write (iout, '(3a)') "24,'", trim(fname_radix)//"_BLZTRP", ".halltens', 'unknown', 'formatted',0" 5207 write (iout, '(3a)') "25,'", trim(fname_radix)//"_BLZTRP", ".trace_fixdoping', 'unknown', 'formatted',0" 5208 write (iout, '(3a)') "26,'", trim(fname_radix)//"_BLZTRP", ".condtens_fixdoping', 'unknown', 'formatted',0" 5209 write (iout, '(3a)') "27,'", trim(fname_radix)//"_BLZTRP", ".halltens_fixdoping', 'unknown', 'formatted',0" 5210 write (iout, '(3a)') "30,'", trim(fname_radix)//"_BLZTRP", "_BZ.dx', 'unknown', 'formatted',0" 5211 write (iout, '(3a)') "31,'", trim(fname_radix)//"_BLZTRP", "_fermi.dx', 'unknown', 'formatted',0" 5212 write (iout, '(3a)') "32,'", trim(fname_radix)//"_BLZTRP", "_sigxx.dx', 'unknown', 'formatted',0" 5213 write (iout, '(3a)') "33,'", trim(fname_radix)//"_BLZTRP", "_sigyy.dx', 'unknown', 'formatted',0" 5214 write (iout, '(3a)') "34,'", trim(fname_radix)//"_BLZTRP", "_sigzz.dx', 'unknown', 'formatted',0" 5215 write (iout, '(3a)') "35,'", trim(fname_radix)//"_BLZTRP", "_band.dat', 'unknown', 'formatted',0" 5216 write (iout, '(3a)') "36,'", trim(fname_radix)//"_BLZTRP", "_band.gpl', 'unknown', 'formatted',0" 5217 write (iout, '(3a)') "37,'", trim(fname_radix)//"_BLZTRP", "_deriv.dat', 'unknown', 'formatted',0" 5218 write (iout, '(3a)') "38,'", trim(fname_radix)//"_BLZTRP", "_mass.dat', 'unknown', 'formatted',0" 5219 5220 close(iout) 5221 end do !isppol 5222 5223 !file is for geometry symmetries etc 5224 filename= trim(fname_radix)//"_BLZTRP.struct" 5225 if (open_file(filename, msg, newunit=iout, form='formatted') /= 0) then 5226 ABI_ERROR(msg) 5227 end if 5228 5229 write (iout, '(a)') "BoltzTraP geometry file generated by ABINIT." 5230 5231 !here we need to print out the unit cell vectors 5232 write (iout, '(3E20.10)') crystal%rprimd(:,1) 5233 write (iout, '(3E20.10)') crystal%rprimd(:,2) 5234 write (iout, '(3E20.10)') crystal%rprimd(:,3) 5235 write (iout, '(I7)') crystal%nsym 5236 5237 do isym=1,crystal%nsym 5238 write (iout,'(3(3I5,2x), a, I5)') & 5239 & crystal%symrel(1,:,isym), & 5240 & crystal%symrel(2,:,isym), & 5241 & crystal%symrel(3,:,isym), & 5242 & ' ! symmetry rotation matrix isym = ', isym 5243 end do 5244 5245 close (iout) 5246 5247 ! second file is for eigenvalues 5248 ! two file names for each spin, if necessary 5249 do isppol=1,nsppol 5250 filename=trim(fname_radix)//"_BLZTRP."//spinsuffix(isppol)//"energy"//trim(so_suffix) 5251 5252 if (open_file (filename, msg, newunit=iout, form='formatted') /= 0) then 5253 ABI_ERROR(msg) 5254 end if 5255 5256 write (iout, '(a,I5)') "BoltzTraP eigen-energies file generated by ABINIT. ispin = ", isppol 5257 write (iout, '(I7, I7, E20.10, a)') & 5258 & ebands%nkpt, nsppol, ha2ryd*ebands%fermie, ' ! nk, nspin, Fermi level(Ry) : energies below in Ry' 5259 5260 do ikpt=1,ebands%nkpt 5261 ! these need to be in reduced coordinates 5262 write (iout, '(3E20.10, I7, a)') & 5263 & ebands%kptns(1,ikpt), ebands%kptns(2,ikpt), ebands%kptns(3,ikpt), nband, ' ! kpt nband' 5264 do iband=1,nband 5265 ! output in Rydberg 5266 write (iout, '(E20.10)') ha2ryd*ebands%eig(iband, ikpt, isppol) 5267 end do 5268 end do 5269 5270 close (iout) 5271 end do 5272 5273 !this file is for tau_k 5274 if (present (tau_k)) then 5275 do isppol = 1, nsppol 5276 filename= trim(fname_radix)//"_"//spinsuffix(isppol)//"BLZTRP.tau_k" 5277 if (open_file(filename, msg, newunit=iout, form='formatted') /= 0) then 5278 ABI_ERROR(msg) 5279 end if 5280 5281 write (iout, '(a)') "BoltzTraP tau_k file generated by ANADDB." 5282 write (iout, '(I7, I7, E20.10, a)')& 5283 & ebands%nkpt, nsppol, ha2ryd*ebands%fermie, ' ! nk, nspin, Fermi level(Ry) : energies below in Ry' 5284 5285 do ikpt=1,ebands%nkpt 5286 ! these need to be in reduced coordinates 5287 write (iout, '(3E20.10, I7, a)') & 5288 & ebands%kptns(1,ikpt), ebands%kptns(2,ikpt), ebands%kptns(3,ikpt), nband, ' ! kpt nband' 5289 do iband=1,nband 5290 ! output in eV 5291 write (iout, '(E20.10)') tau_k(isppol,ikpt,iband) 5292 end do 5293 end do 5294 close (iout) 5295 end do 5296 5297 end if 5298 5299 end subroutine ebands_prtbltztrp
m_ebands/ebands_prtbltztrp_tau_out [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_prtbltztrp_tau_out
FUNCTION
output files for BoltzTraP code, which integrates Boltzmann transport quantities over the Fermi surface for different T and chemical potentials. Abinit provides all necessary input files: struct, energy, input file, and def file for the unit definitions of fortran files in BT. See http://www.icams.de/content/departments/ams/madsen/boltztrap.html Output T-depedent tau_k, modified from ebands_prtbltztrp
INPUTS
eigen(mband*nkpt*nsppol) = array for holding eigenvalues (hartree) fermie = Fermi level fname_radix = radix of file names for output nband = number of bands nkpt = number of k points. nsppol = 1 for unpolarized, 2 for spin-polarized nsym = number of symmetries in space group rprimd(3,3) = dimensional primitive translations for real space (bohr) symrel = symmetry operations in reduced coordinates, real space to be used in future xred(3,natom) = reduced dimensionless atomic coordinates
OUTPUT
(only writing, printing)
SOURCE
5331 subroutine ebands_prtbltztrp_tau_out (eigen, tempermin, temperinc, ntemper, fermie, fname_radix, kpt, & 5332 nband, nelec, nkpt, nspinor, nsppol, nsym, rprimd, symrel, tau_k) 5333 5334 !Arguments ------------------------------------ 5335 !scalars 5336 integer, intent(in) :: nsym, nband, nkpt, nsppol, nspinor, ntemper 5337 real(dp), intent(in) :: tempermin, temperinc 5338 real(dp), intent(in) :: nelec 5339 character(len=fnlen), intent(in) :: fname_radix 5340 !arrays 5341 real(dp), intent(in) :: fermie(ntemper) 5342 integer, intent(in) :: symrel(3,3,nsym) 5343 real(dp), intent(in) :: kpt(3,nkpt) 5344 real(dp), intent(in) :: eigen(nband, nkpt, nsppol) 5345 real(dp), intent(in) :: rprimd(3,3) 5346 real(dp), intent(in) :: tau_k(ntemper,nsppol,nkpt,nband) 5347 5348 !Local variables------------------------------- 5349 !scalars 5350 integer :: iout, isym, iband, isppol, ikpt, itemp 5351 real(dp) :: Temp 5352 real(dp),parameter :: ha2ryd = two 5353 character(len=500) :: msg 5354 character(len=fnlen) :: filename,appendix 5355 5356 ! ************************************************************************* 5357 5358 !input file for boltztrap: general info, Ef, Nelec, etc... 5359 do itemp = 1, ntemper 5360 write(appendix,"(i0)") itemp 5361 filename= trim(fname_radix)//"_BLZTRP.intrans_"//trim(appendix) 5362 if (open_file(filename, msg, newunit=iout, form="formatted", action="write") /= 0) then 5363 ABI_ERROR(msg) 5364 end if 5365 5366 write (iout, '(a)') "GENE # Format of input: generic format, with Symmetries" 5367 write (iout, '(a)') "0 0 0 0.0 # iskip (not presently used) idebug setgap shiftgap" 5368 write (iout, '(E15.5,a,F10.4,a)') fermie(itemp)*two, " 0.0005 0.4 ", nelec, & 5369 " # Fermilevel (Ry), energy grid spacing, energy span around Fermilevel, number of electrons" 5370 write (iout, '(a)') "CALC # CALC (calculate expansion coeff), NOCALC read from file" 5371 write (iout, '(a)') "3 # lpfac, number of latt-points per k-point" 5372 write (iout, '(a)') "BOLTZ # run mode (only BOLTZ is supported)" 5373 write (iout, '(a)') ".15 # (efcut) energy range of chemical potential" 5374 write (iout, '(2f8.2,a)')& 5375 tempermin+temperinc*dble(itemp),tempermin+temperinc*dble(itemp), " # Tmax, temperature grid spacing" 5376 write (iout, '(2a)') "-1 # energyrange of bands given ",& 5377 "individual DOS output sig_xxx and dos_xxx (xxx is band number)" 5378 write (iout, '(a)') "TETRA # DOS calculation method. Other possibility is TETRA" 5379 write (iout, '(a)') "No # not using model for relaxation time" 5380 write (iout, '(a)') "3 # Number of doping levels coefficients will be output for" 5381 write (iout, '(a)') "-1.e16 0.0d0 1.e16 # Values of doping levels (in carriers / cm^3" 5382 close(iout) 5383 end do 5384 5385 !files file, with association of all units for Boltztrap 5386 filename= trim(fname_radix)//"_BLZTRP.def" 5387 if (open_file(filename, msg, newunit=iout, form="formatted", action="write") /= 0) then 5388 ABI_ERROR(msg) 5389 end if 5390 write (iout, '(3a)') "5, '", trim(fname_radix)//"_BLZTRP", ".intrans', 'old', 'formatted',0" 5391 write (iout, '(3a)') "6, '", trim(fname_radix)//"_BLZTRP", ".outputtrans', 'unknown', 'formatted',0" 5392 write (iout, '(3a)') "20,'", trim(fname_radix)//"_BLZTRP", ".struct', 'old', 'formatted',0" 5393 if (nspinor == 1) then 5394 write (iout, '(3a)') "10,'", trim(fname_radix)//"_BLZTRP", ".energy', 'old', 'formatted',0" 5395 else if (nspinor == 2) then 5396 write (iout, '(3a)') "10,'", trim(fname_radix)//"_BLZTRP", ".energyso', 'old', 'formatted',0" 5397 end if 5398 write (iout, '(3a)') "10,'", trim(fname_radix)//"_BLZTRP", ".energy', 'old', 'formatted',0" 5399 write (iout, '(3a)') "11,'", trim(fname_radix)//"_BLZTRP", ".tau_k', 'old', 'formatted',0" 5400 write (iout, '(3a)') "48,'", trim(fname_radix)//"_BLZTRP", ".engre', 'unknown', 'unformatted',0" 5401 write (iout, '(3a)') "49,'", trim(fname_radix)//"_BLZTRP", ".transdos', 'unknown', 'formatted',0" 5402 write (iout, '(3a)') "50,'", trim(fname_radix)//"_BLZTRP", ".sigxx', 'unknown', 'formatted',0" 5403 write (iout, '(3a)') "51,'", trim(fname_radix)//"_BLZTRP", ".sigxxx', 'unknown', 'formatted',0" 5404 write (iout, '(3a)') "21,'", trim(fname_radix)//"_BLZTRP", ".trace', 'unknown', 'formatted',0" 5405 write (iout, '(3a)') "22,'", trim(fname_radix)//"_BLZTRP", ".condtens', 'unknown', 'formatted',0" 5406 write (iout, '(3a)') "24,'", trim(fname_radix)//"_BLZTRP", ".halltens', 'unknown', 'formatted',0" 5407 write (iout, '(3a)') "25,'", trim(fname_radix)//"_BLZTRP", ".trace_fixdoping', 'unknown', 'formatted',0" 5408 write (iout, '(3a)') "26,'", trim(fname_radix)//"_BLZTRP", ".condtens_fixdoping', 'unknown', 'formatted',0" 5409 write (iout, '(3a)') "27,'", trim(fname_radix)//"_BLZTRP", ".halltens_fixdoping', 'unknown', 'formatted',0" 5410 write (iout, '(3a)') "30,'", trim(fname_radix)//"_BLZTRP", "_BZ.dx', 'unknown', 'formatted',0" 5411 write (iout, '(3a)') "31,'", trim(fname_radix)//"_BLZTRP", "_fermi.dx', 'unknown', 'formatted',0" 5412 write (iout, '(3a)') "32,'", trim(fname_radix)//"_BLZTRP", "_sigxx.dx', 'unknown', 'formatted',0" 5413 write (iout, '(3a)') "33,'", trim(fname_radix)//"_BLZTRP", "_sigyy.dx', 'unknown', 'formatted',0" 5414 write (iout, '(3a)') "34,'", trim(fname_radix)//"_BLZTRP", "_sigzz.dx', 'unknown', 'formatted',0" 5415 write (iout, '(3a)') "35,'", trim(fname_radix)//"_BLZTRP", "_band.dat', 'unknown', 'formatted',0" 5416 write (iout, '(3a)') "36,'", trim(fname_radix)//"_BLZTRP", "_band.gpl', 'unknown', 'formatted',0" 5417 write (iout, '(3a)') "37,'", trim(fname_radix)//"_BLZTRP", "_deriv.dat', 'unknown', 'formatted',0" 5418 write (iout, '(3a)') "38,'", trim(fname_radix)//"_BLZTRP", "_mass.dat', 'unknown', 'formatted',0" 5419 close(iout) 5420 5421 !file is for geometry symmetries etc 5422 filename= trim(fname_radix)//"_BLZTRP.struct" 5423 if (open_file(filename, msg, newunit=iout, form="formatted", action="write") /= 0) then 5424 ABI_ERROR(msg) 5425 end if 5426 write (iout, '(a)') "BoltzTraP geometry file generated by ABINIT." 5427 5428 !here we need to print out the unit cell vectors 5429 write (iout, '(3E20.10)') rprimd(:,1) 5430 write (iout, '(3E20.10)') rprimd(:,2) 5431 write (iout, '(3E20.10)') rprimd(:,3) 5432 write (iout, '(I7)') nsym 5433 5434 do isym=1, nsym 5435 write (iout,'(3(3I5,2x), a, I5)') & 5436 symrel(1,:,isym), symrel(2,:,isym), symrel(3,:,isym), ' ! symmetry rotation matrix isym = ', isym 5437 end do 5438 close (iout) 5439 5440 !second file is for eigenvalues 5441 if (nspinor == 1) then 5442 filename= trim(fname_radix)//"_BLZTRP.energy" 5443 else if (nspinor == 2) then 5444 filename= trim(fname_radix)//"_BLZTRP.energyso" 5445 end if 5446 5447 if (open_file(filename, msg, newunit=iout, form="formatted", action="write") /= 0) then 5448 ABI_ERROR(msg) 5449 end if 5450 write (iout, '(a)') "BoltzTraP eigen-energies file generated by ABINIT." 5451 write (iout, '(I7, I7, E20.10, a)') nkpt, nsppol, ha2ryd*fermie(1), ' ! nk, nspin, Fermi level(Ry) : energies below in Ry' 5452 do isppol = 1, nsppol 5453 do ikpt = 1, nkpt 5454 ! these need to be in reduced coordinates 5455 write (iout, '(3E20.10, I7, a)') kpt(1,ikpt), kpt(2,ikpt), kpt(3,ikpt), nband, ' ! kpt nband' 5456 do iband = 1, nband 5457 ! output in eV 5458 write (iout, '(E20.10)') ha2ryd*eigen(iband, ikpt, isppol) 5459 end do 5460 end do 5461 end do 5462 close (iout) 5463 5464 !this file is for tau_k 5465 do itemp = 1, ntemper 5466 Temp=tempermin+temperinc*dble(itemp) 5467 5468 write(appendix,"(i0)") itemp 5469 filename= trim(fname_radix)//"_BLZTRP.tau_k_"//trim(appendix) 5470 if (open_file(filename, msg, newunit=iout, form="formatted", action="write") /= 0) then 5471 ABI_ERROR(msg) 5472 end if 5473 write (iout, '(a,f12.6)') "BoltzTraP tau_k file generated by ANADDB for T= ", Temp 5474 write (iout, '(I7, I7, E20.10, a)') nkpt, nsppol, ha2ryd*fermie(itemp), & 5475 ' ! nk, nspin, Fermi level(Ry) : energies below in Ry' 5476 do isppol = 1, nsppol 5477 do ikpt = 1, nkpt 5478 ! these need to be in reduced coordinates 5479 write (iout, '(3E20.10, I7, a)') kpt(1,ikpt), kpt(2,ikpt), kpt(3,ikpt), nband, ' ! kpt nband' 5480 do iband = 1, nband 5481 ! output in sec 5482 write (iout, '(E20.10)') tau_k(itemp,isppol,ikpt,iband) 5483 end do 5484 end do 5485 end do 5486 close (iout) 5487 end do 5488 5489 end subroutine ebands_prtbltztrp_tau_out
m_ebands/ebands_report_gap [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_report_gap
FUNCTION
Print info on the fundamental and direct gap.
INPUTS
ebands<ebands_t>=Info on the band structure, the smearing technique and the physical temperature used. [header]=Optional title. [unit]=Optional unit for output (std_out if not specified) [mode_paral]=Either "COLL" or "PERS", former is default.
OUTPUT
writing.
[gaps(3,nsppol)]=Fundamental and direct gaps. The third index corresponds to a "status":
0.0dp if gaps were not computed (because there are only valence bands);
-1.0dp if the system (or spin-channel) is metallic;
1.0dp if the gap has been computed.
SOURCE
2800 subroutine ebands_report_gap(ebands, header, unit, mode_paral, gaps) 2801 2802 !Arguments ------------------------------------ 2803 !scalars 2804 integer,intent(in),optional :: unit 2805 character(len=4),intent(in),optional :: mode_paral 2806 character(len=*),intent(in),optional :: header 2807 class(ebands_t),intent(in) :: ebands 2808 !arrays 2809 real(dp),optional,intent(out) :: gaps(3,ebands%nsppol) 2810 2811 !Local variables------------------------------- 2812 !scalars 2813 integer :: ikibz,nband_k,spin,nsppol,ikopt,ivk,ick,ivb,icb,unt,first 2814 real(dp),parameter :: tol_fermi = tol6 2815 real(dp) :: fun_gap,opt_gap 2816 logical :: ismetal 2817 character(len=4) :: my_mode 2818 character(len=500) :: msg 2819 !arrays 2820 integer :: val_idx(ebands%nkpt,ebands%nsppol) 2821 real(dp) :: top_valence(ebands%nkpt),bot_conduct(ebands%nkpt) 2822 2823 ! ********************************************************************* 2824 2825 nsppol = ebands%nsppol 2826 2827 unt =std_out; if (PRESENT(unit )) unt =unit 2828 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral 2829 2830 if (PRESENT(gaps)) gaps=zero 2831 2832 val_idx(:,:) = ebands_get_valence_idx(ebands, tol_fermi) 2833 first = 0 2834 2835 ! Initialize the return status for the gaps 2836 if (PRESENT(gaps)) gaps(1:3,1:nsppol)=zero 2837 2838 do spin=1,nsppol 2839 2840 ! No output if system i metallic 2841 ismetal=ANY(val_idx(:,spin)/=val_idx(1,spin)) 2842 if (ismetal .or. (ebands%occopt==9)) then 2843 if (PRESENT(gaps)) gaps(3,nsppol)=-one 2844 cycle 2845 endif 2846 2847 first=first+1 2848 if (first==1) then 2849 msg=ch10 2850 if (PRESENT(header)) msg=ch10//' === '//TRIM(ADJUSTL(header))//' === ' 2851 call wrtout(unt,msg,my_mode) 2852 end if 2853 2854 ivb=val_idx(1,spin) 2855 icb=ivb+1 2856 2857 do ikibz=1,ebands%nkpt 2858 nband_k = ebands%nband(ikibz+(spin-1)*ebands%nkpt) 2859 top_valence(ikibz) = ebands%eig(ivb,ikibz,spin) 2860 if (icb>nband_k) then 2861 GOTO 10 ! Only occupied states are present, no output! 2862 end if 2863 bot_conduct(ikibz) = ebands%eig(icb,ikibz,spin) 2864 end do 2865 2866 ! Get minimum of the direct Gap 2867 ikopt= imin_loc(bot_conduct-top_valence) 2868 opt_gap=bot_conduct(ikopt)-top_valence(ikopt) 2869 2870 ! Get fundamental Gap 2871 ick = imin_loc(bot_conduct) 2872 ivk = imax_loc(top_valence) 2873 fun_gap = ebands%eig(icb,ick,spin)-ebands%eig(ivb,ivk,spin) 2874 2875 write(msg,'(a,i2,a,2(a,f8.4,a,3f8.4,a),33x,a,3f8.4)')& 2876 ' >>>> For spin ',spin,ch10,& 2877 ' Minimum direct gap = ',opt_gap*Ha_eV,' [eV], located at k-point : ',ebands%kptns(:,ikopt),ch10,& 2878 ' Fundamental gap = ',fun_gap*Ha_eV,' [eV], Top of valence bands at : ',ebands%kptns(:,ivk),ch10, & 2879 ' Bottom of conduction at : ',ebands%kptns(:,ick) 2880 call wrtout(unt,msg,my_mode) 2881 2882 if (present(gaps)) gaps(:,spin) = [fun_gap, opt_gap, one] 2883 end do !spin 2884 2885 return 2886 2887 10 continue 2888 call wrtout(std_out, "Not enough states to calculate the band gap.", "COLL") 2889 2890 end subroutine ebands_report_gap
m_ebands/ebands_set_extrael [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_set_extrael
FUNCTION
Add extrael to ebands%nelect. Set value of ebands%extrael Recompute Fermi level from eigenenergies and new occupation numbers according to the smearing scheme defined by occopt (and smearing width tsmear or tphysel) as well as entropy and derivative of occupancies wrt the energy for each band and k point.
INPUTS
nelect=New number of electrons nholes=New number of excited holes extrael=Number of electrons be added in units. Negative to add holes spinmagntarget=if differ from -99.99d0, fix the spin polarization (in Bohr magneton) [prtvol]=Verbosity level
OUTPUT
msg=String describing the changes in fermie and nelect.
NOTES
The routine assumes metallic occupation scheme and will abort it this condition is not satisfied. Use ebands_set_scheme before calling this routine, if you have a semiconductor.
SOURCE
2595 subroutine ebands_set_extrael(ebands, nelect, nholes, spinmagntarget, msg, prtvol) 2596 2597 !Arguments ------------------------------------ 2598 !scalars 2599 class(ebands_t),intent(inout) :: ebands 2600 integer,optional,intent(in) :: prtvol 2601 real(dp),intent(in) :: nelect,nholes,spinmagntarget 2602 character(len=*),intent(out) :: msg 2603 2604 !Local variables------------------------------- 2605 !scalars 2606 integer :: my_prtvol 2607 real(dp) :: prev_fermie,prev_fermih,prev_nelect,prev_nholes 2608 2609 ! ************************************************************************* 2610 2611 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 2612 2613 if (.not. ebands_has_metal_scheme(ebands)) then 2614 ABI_ERROR("set_extrael assumes a metallic occupation scheme. Use ebands_set_scheme!") 2615 end if 2616 2617 prev_fermie = ebands%fermie; prev_nelect = ebands%nelect 2618 prev_fermih = ebands%fermie; prev_nholes = zero 2619 ! Here we set the value of extrael 2620 ebands%extrael = nelect-nholes 2621 ebands%nelect = ebands%nelect + ebands%extrael 2622 if (ebands%occopt /= 9) then 2623 ebands%ne_qFD = zero 2624 ebands%nh_qFD = zero 2625 else 2626 prev_fermie = ebands%fermie; prev_nelect = ebands%ne_qFD 2627 prev_fermih = ebands%fermih; prev_nholes = ebands%nh_qFD 2628 ebands%ne_qFD = nelect 2629 ebands%nh_qFD = nholes 2630 end if 2631 2632 call ebands_update_occ(ebands, spinmagntarget, prtvol=my_prtvol) 2633 2634 if (ebands%occopt/=9) then 2635 write(msg,"(2(a,es16.6),a,2(a,es16.6))")& 2636 " Old fermi level: ",prev_fermie,", with nelect: ",prev_nelect,ch10,& 2637 " New fermi level: ",ebands%fermie,", with nelect: ",ebands%nelect 2638 call wrtout(std_out, msg) 2639 else 2640 write(msg,"(2(a,es16.6),a,2(a,es16.6))")& 2641 " Old electron fermi level: ",prev_fermie,", with nelect: ",prev_nelect,ch10,& 2642 " New electron fermi level: ",ebands%fermie,", with nelect: ",ebands%ne_qFD 2643 call wrtout(std_out, msg) 2644 write(msg,"(2(a,es16.6),a,2(a,es16.6))")& 2645 " Old holes fermi level: ",prev_fermih,", with nelect: ",prev_nelect-prev_nholes,ch10,& 2646 " New holes fermi level: ",ebands%fermih,", with nholes: ",ebands%nelect-ebands%nh_qFD 2647 call wrtout(std_out, msg) 2648 end if 2649 2650 end subroutine ebands_set_extrael
m_ebands/ebands_set_fermie [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_set_fermie
FUNCTION
Set the new Fermi level from eigenenergies eigen and change the number of electrons Compute also new occupation numbers at each k point, from eigenenergies eigen, according to the smearing scheme defined by occopt (and smearing width tsmear or tphysel) as well as entropy and derivative of occupancies wrt the energy for each band and k point.
INPUTS
fermie=New fermi level
OUTPUT
msg=String describing the changes in fermie and nelect.
NOTES
The routine assumes metallic occupation scheme and will abort it this condition is not satisfied. Use ebands_set_scheme before calling this routine, if you have a semiconductor.
SOURCE
2507 subroutine ebands_set_fermie(ebands, fermie, msg) 2508 2509 !Arguments ------------------------------------ 2510 !scalars 2511 class(ebands_t),intent(inout) :: ebands 2512 real(dp),intent(in) :: fermie 2513 character(len=*),intent(out) :: msg 2514 2515 !Local variables------------------------------- 2516 !scalars 2517 integer,parameter :: option1=1, unitdos0 = 0 2518 integer :: mband, nkpt, nsppol 2519 real(dp),parameter :: dosdeltae0 = zero 2520 real(dp) :: prev_fermie, prev_nelect, maxocc 2521 !arrays 2522 real(dp),allocatable :: doccde(:),occ(:),eigen(:) 2523 2524 ! ************************************************************************* 2525 2526 if (ebands%occopt == 9) then 2527 ABI_ERROR("set_fermie unavailable when occopt 9") 2528 end if 2529 if (.not. ebands_has_metal_scheme(ebands)) then 2530 ABI_ERROR("set_fermie assumes a metallic occupation scheme. Use ebands_set_scheme before calling ebands_set_fermie!") 2531 end if 2532 2533 prev_fermie = ebands%fermie; prev_nelect = ebands%nelect 2534 2535 ! newocc assumes eigenvalues and occupations packed in 1d-vector!! 2536 mband = ebands%mband 2537 nkpt = ebands%nkpt 2538 nsppol = ebands%nsppol 2539 maxocc = two / (nsppol*ebands%nspinor) 2540 2541 ABI_MALLOC(eigen, (mband*nkpt*nsppol)) 2542 call get_eneocc_vect(ebands, 'eig', eigen) 2543 ABI_MALLOC(occ, (mband*nkpt*nsppol)) 2544 ABI_MALLOC(doccde, (mband*nkpt*nsppol)) 2545 2546 ! Get the total number of electrons nelect, given the new fermi energy. 2547 call getnel(doccde,dosdeltae0,eigen,ebands%entropy,fermie,fermie,maxocc,mband,ebands%nband,& 2548 ebands%nelect,nkpt,nsppol,occ,ebands%occopt,option1,ebands%tphysel,ebands%tsmear,unitdos0,& 2549 ebands%wtk,1,ebands%nband(1)) 2550 2551 ! Save changes in ebands%. 2552 ebands%fermie = fermie 2553 call put_eneocc_vect(ebands,'occ' ,occ) 2554 call put_eneocc_vect(ebands,'doccde',doccde) 2555 2556 ABI_FREE(eigen) 2557 ABI_FREE(occ) 2558 ABI_FREE(doccde) 2559 2560 write(msg,"(2(a,es16.6),a,2(a,es16.6))") & 2561 " Old fermi level: ",prev_fermie,", with nelect: ",prev_nelect,ch10, & 2562 " New fermi level: ",ebands%fermie,", with nelect: ",ebands%nelect 2563 2564 end subroutine ebands_set_fermie
m_ebands/ebands_set_scheme [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_set_scheme
FUNCTION
Set the occupation scheme and re-calculate new occupation numbers, the Fermi level and the Max occupied band index for each spin channel starting from the the knowledge of eigenvalues. See ebands_update_occ for more info.
INPUTS
occopt=Occupation options (see input variable) tsmear=Temperature of smearing. spinmagntarget=if differ from -99.99d0, fix the spin polarization (in Bohr magneton) prtvol=Verbosity level (0 for lowest level) [update_occ]=False to avoid recomputing occupation factors (mainly used when a call to set_scheme is followed by another call to set_extrael (update_occ is expensive for large k-meshes). Default: True.
SOURCE
2449 subroutine ebands_set_scheme(ebands, occopt, tsmear, spinmagntarget, prtvol, update_occ) 2450 2451 !Arguments ------------------------------------ 2452 !scalars 2453 class(ebands_t),intent(inout) :: ebands 2454 integer,intent(in) :: occopt 2455 integer,intent(in) :: prtvol 2456 real(dp),intent(in) :: tsmear, spinmagntarget 2457 logical,optional,intent(in) :: update_occ 2458 2459 !Local variables------------------------------- 2460 !scalars 2461 logical :: my_update_occ 2462 2463 ! ************************************************************************* 2464 2465 my_update_occ = .True.; if (present(update_occ)) my_update_occ = update_occ 2466 2467 if (prtvol > 10) then 2468 call wrtout(std_out, " Changing occupation scheme in electron bands") 2469 call wrtout(std_out, sjoin(" occopt:", itoa(ebands%occopt), " ==> ", itoa(occopt))) 2470 call wrtout(std_out, sjoin(" tsmear:", ftoa(ebands%tsmear), " ==> ", ftoa(tsmear))) 2471 end if 2472 2473 ebands%occopt = occopt; ebands%tsmear = tsmear 2474 2475 if (my_update_occ) then 2476 call ebands_update_occ(ebands, spinmagntarget, prtvol=prtvol) 2477 if (prtvol > 10) call wrtout(std_out, sjoin(' Fermi level is now:', ftoa(ebands%fermie))) 2478 end if 2479 2480 end subroutine ebands_set_scheme
m_ebands/ebands_sort [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_sort
FUNCTION
Sort eigvalues_k in ascending order and reorder arrays depending on nband_k Mainly used when interpolating band energies as the interpolator may not produce ordered eigenvalues and there are routines whose implementation assumes eig(b) <= eig(b+1)
SIDE EFFECTS
ebands<ebands_t> = Object with input energies sorted in output.
SOURCE
4140 subroutine ebands_sort(self) 4141 4142 !Arguments ------------------------------------ 4143 !scalars 4144 class(ebands_t),intent(inout) :: self 4145 4146 !Local variables------------------------------- 4147 !scalars 4148 integer :: spin, ik_ibz, band, nband_k 4149 !arrays 4150 integer :: iperm_k(self%mband) 4151 4152 ! ********************************************************************* 4153 4154 do spin=1,self%nsppol 4155 do ik_ibz=1,self%nkpt 4156 nband_k = self%nband(ik_ibz + (spin - 1) * self%nkpt) 4157 iperm_k = [(band, band=1, nband_k)] 4158 call sort_dp(nband_k, self%eig(:, ik_ibz, spin), iperm_k, tol12) 4159 4160 ! Shuffle other arrays depending on nband_k 4161 self%occ(1:nband_k, ik_ibz, spin) = self%occ(iperm_k(1:nband_k), ik_ibz, spin) 4162 self%doccde(1:nband_k, ik_ibz, spin) = self%doccde(iperm_k(1:nband_k), ik_ibz, spin) 4163 !if (allocated(self%velocity)) then 4164 ! self%velocity(:, 1:nband_k, ik_ibz, spin) = self%velocity(:, iperm_k(1:nband_k), ik_ibz, spin) 4165 !end if 4166 end do 4167 end do 4168 4169 end subroutine ebands_sort
m_ebands/ebands_update_occ [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_update_occ
FUNCTION
Calculate new occupation numbers, the Fermi level and the Max occupied band index for each spin channel starting from the the knowledge of eigenvalues.
INPUTS
ebands<ebands_t>=Info on the band structure, the smearing technique and the physical temperature used. spinmagntarget=if differ from -99.99d0, fix the spin polarization (in Bohr magneton) [stmbias]= [prtvol]=Verbosity level (0 for lowest level) [fermie_to_zero]=If True, fermie is set to zero and energies are shifted accordingly. Default: False
OUTPUT
see also SIDE EFFECTS.
SIDE EFFECTS
=== For metallic occupation the following quantites are recalculated === %fermie=the new Fermi energy %entropy=the new entropy associated with the smearing. %occ(mband,nkpt,nsppol)=occupation numbers %doccde(mband,nkpt,nsppol)=derivative of occupancies wrt the energy for each band and k point === In case of semiconductors === All the quantitities in ebands are left unchanged with the exception of: %fermie=Redefined so that it is in the middle of the gap %entropy=Set to zero
SOURCE
2278 subroutine ebands_update_occ(ebands, spinmagntarget, stmbias, prtvol, fermie_to_zero) 2279 2280 !Arguments ------------------------------------ 2281 !scalars 2282 class(ebands_t),intent(inout) :: ebands 2283 integer,optional,intent(in) :: prtvol 2284 real(dp),intent(in) :: spinmagntarget 2285 real(dp),optional,intent(in) :: stmbias 2286 logical,optional,intent(in) :: fermie_to_zero 2287 2288 !Local variables------------------------------- 2289 !scalars 2290 integer :: band,mband,ikibz,nkpt,spin,nsppol,my_prtvol,nband_k 2291 real(dp) :: entropy,fermie,fermih,stmbias_local,ndiff,cbot,vtop,maxocc 2292 character(len=500) :: msg 2293 !arrays 2294 real(dp) :: nelect_spin(ebands%nsppol),condbottom(ebands%nsppol),valencetop(ebands%nsppol) 2295 real(dp),allocatable :: doccde(:),occ(:),eigen(:) 2296 2297 ! ************************************************************************* 2298 2299 my_prtvol = 0; if (PRESENT(prtvol )) my_prtvol = prtvol 2300 stmbias_local = zero; if (PRESENT(stmbias)) stmbias_local = stmbias 2301 2302 if (ebands_has_metal_scheme(ebands)) then 2303 ! Compute new occupation numbers if metallic occupation. 2304 if (my_prtvol > 10) then 2305 call wrtout(std_out, sjoin(' metallic scheme: calling newocc with spinmagntarget:', ftoa(spinmagntarget, fmt="f9.5"))) 2306 end if 2307 2308 ! newocc assumes eigenvalues and occupations packed in 1d-vector!! 2309 mband = ebands%mband; nkpt = ebands%nkpt; nsppol = ebands%nsppol 2310 2311 ABI_MALLOC(eigen, (mband*nkpt*nsppol)) 2312 call get_eneocc_vect(ebands, 'eig', eigen) 2313 2314 ABI_MALLOC(occ, (mband*nkpt*nsppol)) 2315 ABI_MALLOC(doccde, (mband*nkpt*nsppol)) 2316 2317 call newocc(doccde,eigen,entropy,fermie,fermih,ebands%ivalence,spinmagntarget,mband,ebands%nband,& 2318 ebands%nelect,ebands%ne_qFD,ebands%nh_qFD,ebands%nkpt,ebands%nspinor,ebands%nsppol,occ,ebands%occopt,& 2319 my_prtvol,ebands%tphysel,ebands%tsmear,ebands%wtk,stmbias=stmbias_local) 2320 2321 ! Save output in ebands%. 2322 ebands%entropy = entropy 2323 ebands%fermie = fermie 2324 ebands%fermih = fermih 2325 call put_eneocc_vect(ebands, 'occ', occ) 2326 call put_eneocc_vect(ebands, 'doccde', doccde) 2327 ABI_FREE(eigen) 2328 ABI_FREE(occ) 2329 ABI_FREE(doccde) 2330 2331 else 2332 ! Semiconductor (non magnetic case) 2333 maxocc = two / (ebands%nsppol*ebands%nspinor) 2334 ! 2335 ! FIXME here there is an inconsistency btw the GW code and Abinit 2336 ! In ABINIT, Fermi is set to the HOMO level while in GW fermi is at midgap 2337 ! In case of crystal systems, the later convention should be preferable. 2338 ! Anyway we have to decide and follow a unique convention to avoid problems. 2339 ! 2340 ! Occupation factors MUST be initialized 2341 !if (ALL(ABS(ebands%occ) < tol6)) then 2342 if (ebands%occopt /= 2) then 2343 ABI_CHECK(ebands%nelect == nint(ebands%nelect), "nelect should be integer") 2344 mband = nint((ebands%nelect * ebands%nspinor) / 2) 2345 ebands%occ = zero 2346 ebands%occ(1:mband,:,:) = maxocc 2347 !ABI_ERROR("Occupation factors are not initialized, likely due to scf = -2") 2348 end if 2349 2350 ! Calculate the valence index for each spin channel. 2351 do spin=1,ebands%nsppol 2352 valencetop(spin) = smallest_real 2353 condbottom(spin) = greatest_real 2354 do ikibz=1,ebands%nkpt 2355 nband_k = ebands%nband(ikibz + (spin-1)*ebands%nkpt) 2356 do band=1,nband_k 2357 if (ebands%occ(band,ikibz,spin) / maxocc > one-tol6 .and. valencetop(spin) < ebands%eig(band,ikibz,spin)) then 2358 valencetop(spin) = ebands%eig(band,ikibz,spin) 2359 end if 2360 if (ebands%occ(band,ikibz,spin) / maxocc < tol6 .and. condbottom(spin) > ebands%eig(band,ikibz,spin)) then 2361 condbottom(spin) = ebands%eig(band,ikibz,spin) 2362 end if 2363 end do 2364 end do 2365 end do 2366 2367 vtop = MAXVAL(valencetop) 2368 cbot = MINVAL(condbottom) 2369 2370 write(msg,'(3(a,f8.4,2a))') & 2371 " Top of valence: ", vtop * Ha_eV," (eV)", ch10, & 2372 " Bottom of conduction: ", cbot * Ha_eV, " (eV)", ch10, & 2373 " Fundamental gap:", (cbot - vtop) * Ha_eV, " (eV)", ch10 2374 call wrtout(std_out, msg) 2375 2376 if (ebands%nsppol == 2) then 2377 if (ABS(vtop - MINVAL(valencetop)) > tol6) then 2378 call wrtout(std_out, sjoin(' Top of valence is spin: ', itoa(imax_loc(valencetop)))) 2379 end if 2380 if (ABS(cbot - MAXVAL(condbottom)) > tol6) then 2381 call wrtout(std_out, ' Bottom of conduction is spin: ', itoa(imin_loc(condbottom))) 2382 end if 2383 end if 2384 2385 ! Save results. Here I dont know if it is better to be consistent with the abinit convention i.e fermi = vtop 2386 ebands%entropy = zero 2387 ebands%fermie = (vtop + cbot) / 2 2388 if (ABS(cbot - vtop) < tol4) ebands%fermie = vtop ! To avoid error on the last digit 2389 end if 2390 2391 call wrtout(std_out, sjoin(' Fermi level: ', ftoa(ebands%fermie * Ha_eV, fmt="f8.4"), " (eV)", ch10)) 2392 2393 ! Compute number of electrons for each spin channel. 2394 nelect_spin(:)=zero 2395 do spin=1,ebands%nsppol 2396 do ikibz=1,ebands%nkpt 2397 nband_k = ebands%nband(ikibz+(spin-1)*ebands%nkpt) 2398 nelect_spin(spin)= nelect_spin(spin) + ebands%wtk(ikibz) * sum(ebands%occ(1:nband_k,ikibz,spin)) 2399 end do 2400 end do 2401 2402 ndiff = ebands%nelect - SUM(nelect_spin) 2403 if (my_prtvol > 0) then 2404 write(msg,'(2a,f6.2,2a,f7.4)')ch10,& 2405 ' Total number of electrons: ', sum(nelect_spin),ch10,& 2406 ' Input and calculated no. of electrons differ by ',ndiff 2407 call wrtout(std_out, msg) 2408 end if 2409 2410 if (ABS(ndiff) > 5.d-2*ebands%nelect) then 2411 write(msg,'(2a,2(a,es12.4))') & 2412 'Too large difference in number of electrons:,',ch10,& 2413 'Expected = ',ebands%nelect,' Calculated = ',sum(nelect_spin) 2414 ABI_ERROR(msg) 2415 end if 2416 2417 if (present(fermie_to_zero)) then 2418 if (fermie_to_zero) then 2419 ebands%eig = ebands%eig - ebands%fermie 2420 ebands%fermih = ebands%fermih - ebands%fermie 2421 ebands%fermie = zero 2422 end if 2423 end if 2424 2425 end subroutine ebands_update_occ
m_ebands/ebands_vcbm_range_from_gaps [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_vcbm_range_from_gaps
FUNCTION
Find band and energy range for states close to the CBM/VBM given input energies in ebands and gaps. Return exit status and error message in msg.
INPUTS
gaps<gaps_t>=Object with info on the gaps. erange(2)=Energy range for holes and electrons. Only those states whose relative position wrt to the VBM/CBM is <= than erange are icluded. Note that relative positions are always positive (even for holes). Use a negative value to exclude either holes or electrons.
OUTPUT
e_lowhigh(2)=min and Max energy. band_lowhigh=min and Max band index. [ks_range]: For each spin and k-point, the min and max band index included in the output set. if (ik, spin) is not included then ib_work(1, ik, spin) > ib_work(2, ik, spin) = -huge(1)
SOURCE
1630 integer function ebands_vcbm_range_from_gaps(ebands, gaps, erange, e_lowhigh, band_lowhigh, ks_range, msg) result(ierr) 1631 1632 !Arguments ------------------------------------ 1633 !scalars 1634 class(ebands_t),intent(in) :: ebands 1635 class(gaps_t),intent(in) :: gaps 1636 real(dp),intent(in) :: erange(2) 1637 real(dp),intent(out) :: e_lowhigh(2) 1638 integer,intent(out) :: band_lowhigh(2) 1639 integer,optional,intent(out) :: ks_range(2, ebands%nkpt, ebands%nsppol) 1640 character(len=*),intent(out) :: msg 1641 1642 !Local variables------------------------------- 1643 integer :: band, ik, spin, band_low, band_high 1644 real(dp) :: cmin, vmax, ee, elow, ehigh 1645 integer,allocatable :: ib_work(:,:,:) 1646 1647 ! ************************************************************************* 1648 1649 ABI_MALLOC(ib_work, (2, ebands%nkpt, ebands%nsppol)) 1650 elow = huge(one); ehigh = -huge(one) 1651 band_low = huge(1); band_high = -huge(1) 1652 1653 ierr = 1 1654 do spin=1,ebands%nsppol 1655 ! Get cmb and vbm with some tolerance 1656 vmax = gaps%vb_max(spin) + tol2 * eV_Ha 1657 cmin = gaps%cb_min(spin) - tol2 * eV_Ha 1658 do ik=1,ebands%nkpt 1659 ib_work(1, ik, spin) = huge(1) 1660 ib_work(2, ik, spin) = -huge(1) 1661 do band=1,ebands%nband(ik+(spin-1)*ebands%nkpt) 1662 ee = ebands%eig(band, ik, spin) 1663 if (erange(1) > zero) then 1664 if (ee <= vmax .and. vmax - ee <= erange(1)) then 1665 ib_work(1, ik, spin) = min(ib_work(1, ik, spin), band) 1666 ib_work(2, ik, spin) = max(ib_work(2, ik, spin), band) 1667 elow = min(elow, ee); ehigh = max(ehigh, ee) 1668 band_low = min(band_low, band); band_high = max(band_high, band) 1669 !write(std_out, *), "Adding valence", band 1670 end if 1671 end if 1672 if (erange(2) > zero) then 1673 if (ee >= cmin .and. ee - cmin <= erange(2)) then 1674 ib_work(1, ik, spin) = min(ib_work(1, ik, spin), band) 1675 ib_work(2, ik, spin) = max(ib_work(2, ik, spin), band) 1676 elow = min(elow, ee); ehigh = max(ehigh, ee) 1677 band_low = min(band_low, band); band_high = max(band_high, band) 1678 !write(std_out, *)"Adding conduction", band 1679 end if 1680 end if 1681 end do 1682 end do 1683 end do 1684 1685 e_lowhigh = [elow, ehigh] 1686 band_lowhigh = [band_low, band_high] 1687 1688 if (present(ks_range)) ks_range = ib_work 1689 ABI_FREE(ib_work) 1690 1691 ! Set exit status and msg. Caller will handle it. 1692 ierr = 0; msg = "" 1693 if (elow > ehigh) then 1694 ierr = 1 1695 write(msg, *)"Cannot find states close to the band edges with erange: ", erange 1696 end if 1697 1698 end function ebands_vcbm_range_from_gaps
m_ebands/ebands_write [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_write
FUNCTION
Driver routine to write bands in different (txt) formats. This routine should be called by a single processor.
INPUTS
prtebands=Flag selecting the output format:
0 --> None
1 --> xmgrace
2 --> gnuplot (not coded yet)
3 --> EIG format (not coded yet)
prefix=Prefix for output filename.
[kptbounds(:,:)]=Optional argument giving the extrema of the k-path.
OUTPUT
Only writing.
SOURCE
5516 subroutine ebands_write(ebands, prtebands, prefix, kptbounds) 5517 5518 !Arguments ------------------------------------ 5519 !scalars 5520 integer,intent(in) :: prtebands 5521 type(ebands_t),intent(in) :: ebands 5522 character(len=*),intent(in) :: prefix 5523 real(dp),optional,intent(in) :: kptbounds(:,:) 5524 5525 ! ********************************************************************* 5526 5527 select case (prtebands) 5528 case (0) 5529 return 5530 case (1) 5531 !call wrtout(std_out, sjoin(" Writing interpolated bands to:", path) 5532 if (present(kptbounds)) then 5533 call ebands_write_xmgrace(ebands, strcat(prefix, "_EBANDS.agr"), kptbounds=kptbounds) 5534 else 5535 call ebands_write_xmgrace(ebands, strcat(prefix, "_EBANDS.agr")) 5536 end if 5537 case (2) 5538 !call wrtout(std_out, sjoin(" Writing interpolated bands to:", path) 5539 if (present(kptbounds)) then 5540 call ebands_write_gnuplot(ebands, prefix, kptbounds=kptbounds) 5541 else 5542 call ebands_write_gnuplot(ebands, prefix) 5543 end if 5544 case default 5545 ABI_WARNING(sjoin("Unsupported value for prtebands:", itoa(prtebands))) 5546 end select 5547 5548 end subroutine ebands_write
m_ebands/ebands_write_bxsf [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_write_bxsf
FUNCTION
Write 3D energies for Fermi surface visualization (XSF format)
INPUTS
ebands<ebands_t>=The object describing the band structure. crystal<crystal_t>=Info on unit cell and symmetries. fname=File name for output.
OUTPUT
ierr=Status error.
SIDE EFFECTS
Produce BXSF file.
SOURCE
2212 integer function ebands_write_bxsf(ebands, crystal, fname) result(ierr) 2213 2214 !Arguments ------------------------------------ 2215 !scalars 2216 character(len=*),intent(in) :: fname 2217 class(ebands_t),intent(in) :: ebands 2218 class(crystal_t),intent(in) :: crystal 2219 2220 !Local variables------------------------------- 2221 logical :: use_timrev 2222 2223 ! ************************************************************************* 2224 2225 use_timrev = (crystal%timrev==2) 2226 2227 if (ebands%occopt /= 9) then 2228 call printbxsf(ebands%eig,zero,ebands%fermie,crystal%gprimd,ebands%kptrlatt,ebands%mband,& 2229 ebands%nkpt,ebands%kptns,crystal%nsym,crystal%use_antiferro,crystal%symrec,crystal%symafm,& 2230 use_timrev,ebands%nsppol,ebands%shiftk,ebands%nshiftk,fname,ierr) 2231 else 2232 call printbxsf(ebands%eig,zero,ebands%fermie,crystal%gprimd,ebands%kptrlatt,ebands%mband,& 2233 ebands%nkpt,ebands%kptns,crystal%nsym,crystal%use_antiferro,crystal%symrec,crystal%symafm,& 2234 use_timrev,ebands%nsppol,ebands%shiftk,ebands%nshiftk,trim(fname)//"-e",ierr) 2235 2236 call printbxsf(ebands%eig,zero,ebands%fermih,crystal%gprimd,ebands%kptrlatt,ebands%mband,& 2237 ebands%nkpt,ebands%kptns,crystal%nsym,crystal%use_antiferro,crystal%symrec,crystal%symafm,& 2238 use_timrev,ebands%nsppol,ebands%shiftk,ebands%nshiftk,trim(fname)//"-h",ierr) 2239 end if 2240 2241 end function ebands_write_bxsf
m_ebands/ebands_write_gnuplot [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_write_gnuplot
FUNCTION
Write bands in gnuplot format. This routine should be called by a single processor. Use the driver `ebands_write` to support different formats.
INPUTS
prefix=prefix for files (.data, .gnuplot) [kptbounds(:,:)]=Optional argument giving the extrema of the k-path.
OUTPUT
Only writing
SOURCE
5695 subroutine ebands_write_gnuplot(ebands, prefix, kptbounds) 5696 5697 !Arguments ------------------------------------ 5698 !scalars 5699 type(ebands_t),intent(in) :: ebands 5700 character(len=*),intent(in) :: prefix 5701 real(dp),optional,intent(in) :: kptbounds(:,:) 5702 5703 !Local variables------------------------------- 5704 !scalars 5705 integer :: unt,gpl_unt,ik,spin,ii,start,nkbounds 5706 character(len=500) :: msg,fmt 5707 character(len=fnlen) :: datafile,basefile 5708 !arrays 5709 integer :: g0(3) 5710 integer,allocatable :: bounds2kpt(:) 5711 5712 ! ********************************************************************* 5713 5714 nkbounds = 0 5715 if (present(kptbounds)) then 5716 if (product(shape(kptbounds)) > 0 ) then 5717 ! Find correspondence between kptbounds and k-points in ebands. 5718 nkbounds = size(kptbounds, dim=2) 5719 ABI_MALLOC(bounds2kpt, (nkbounds)) 5720 bounds2kpt = 1; start = 1 5721 do ii=1,nkbounds 5722 do ik=start,ebands%nkpt 5723 if (isamek(ebands%kptns(:, ik), kptbounds(:, ii), g0)) then 5724 bounds2kpt(ii) = ik; start = ik + 1; exit 5725 end if 5726 end do 5727 end do 5728 end if 5729 end if 5730 5731 datafile = strcat(prefix, "_EBANDS.data") 5732 if (open_file(datafile, msg, newunit=unt, form="formatted", action="write") /= 0) then 5733 ABI_ERROR(msg) 5734 end if 5735 if (open_file(strcat(prefix, "_EBANDS.gnuplot"), msg, newunit=gpl_unt, form="formatted", action="write") /= 0) then 5736 ABI_ERROR(msg) 5737 end if 5738 basefile = basename(datafile) 5739 5740 write(unt,'(a)') "# Electron band structure data file" 5741 write(unt,'(a)') "# Generated by Abinit" 5742 write(unt,'(4(a,i0))') & 5743 "# mband: ",ebands%mband,", nkpt: ",ebands%nkpt,", nsppol: ",ebands%nsppol,", nspinor: ",ebands%nspinor 5744 write(unt,'(a,f8.2,a,i0,2(a,f8.2))') & 5745 "# nelect: ",ebands%nelect,", occopt: ",ebands%occopt,", tsmear: ",ebands%tsmear,", tphysel: ",ebands%tphysel 5746 write(unt,'(a,f8.2,a)') "# Energies are in eV. Zero set to efermi, Previously it was at: ",ebands%fermie * Ha_eV, " [eV]" 5747 write(unt,'(a)')"# List of k-points and their index (C notation i.e. count from 0)" 5748 do ik=1,ebands%nkpt 5749 write(unt, "(a)")sjoin("#", itoa(ik-1), ktoa(ebands%kptns(:,ik))) 5750 end do 5751 5752 fmt = sjoin("(i0,1x,", itoa(ebands%mband), "(es16.8,1x))") 5753 write(unt,'(a)') ' ' 5754 do spin=1,ebands%nsppol 5755 write(unt,'(a,i0)') '# [kpt-index, band_1, band_2 ...] for spin: ',spin 5756 do ik=1,ebands%nkpt 5757 write(unt,fmt) ik-1, (ebands%eig(:, ik, spin) - ebands%fermie) * Ha_eV 5758 end do 5759 write(unt,'(a)') ' ' 5760 end do 5761 5762 ! gnuplot script file 5763 write(gpl_unt,'(a)') '# File to plot phonon bandstructure with gnuplot' 5764 write(gpl_unt,'(a)') "#set terminal postscript eps enhanced color font 'Times-Roman,26' lw 2" 5765 write(gpl_unt,'(a)') '#use the next lines to make a nice figure for a paper' 5766 write(gpl_unt,'(a)') '#set term postscript enhanced eps color lw 0.5 dl 0.5' 5767 write(gpl_unt,'(a)') '#set pointsize 0.275' 5768 write(gpl_unt,'(a)') 'set palette defined ( 0 "blue", 3 "green", 6 "yellow", 10 "red" )' 5769 write(gpl_unt,'(a)') 'unset key' 5770 write(gpl_unt,'(a)') '# can make pointsize smaller (~0.5). Too small and nothing is printed' 5771 write(gpl_unt,'(a)') 'set pointsize 0.8' 5772 write(gpl_unt,'(a)') 'set view 0,0' 5773 write(gpl_unt,'(a,i0,a)') 'set xrange [0:',ebands%nkpt-1,']' 5774 write(gpl_unt,'(2(a,es16.8),a)')& 5775 'set yrange [',minval((ebands%eig - ebands%fermie) * Ha_eV),':',maxval((ebands%eig - ebands%fermie) * Ha_eV),']' 5776 write(gpl_unt,'(a)') 'set xlabel "Momentum"' 5777 write(gpl_unt,'(a)') 'set ylabel "Energy [eV]"' 5778 write(gpl_unt,'(a)') strcat('set title "', replace(basefile, "_", "\\_"), '"') 5779 if (nkbounds == 0) then 5780 write(gpl_unt,'(a)') 'set grid xtics' 5781 else 5782 write(gpl_unt,"(a)")"# Add vertical lines in correspondence of high-symmetry points." 5783 write(gpl_unt,'(a)') 'unset xtics' 5784 do ii=1,nkbounds 5785 write(gpl_unt,"(a,2(i0,a))") & 5786 "set arrow from ",bounds2kpt(ii)-1,",graph(0,0) to ",bounds2kpt(ii)-1,",graph(1,1) nohead ls 'dashed'" 5787 !write(gpl_unt,"(a)")sjoin("set xtics add('kname'", itoa(bounds2kpt(ii)-1), ")") 5788 end do 5789 5790 end if 5791 write(gpl_unt,"(a)")sjoin("mband =", itoa(ebands%mband)) 5792 write(gpl_unt,"(a)")strcat('plot for [i=2:mband] "', basefile, '" u 1:i every :1 with lines linetype -1') 5793 if (ebands%nsppol == 2) then 5794 write(gpl_unt,"(a)")strcat('replot for [i=2:mband] "', basefile, '" u 1:i every :2 with lines linetype 4') 5795 end if 5796 write(gpl_unt, "(a)")"pause -1" 5797 5798 close(unt) 5799 close(gpl_unt) 5800 5801 ABI_SFREE(bounds2kpt) 5802 5803 end subroutine ebands_write_gnuplot
m_ebands/ebands_write_nesting [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_write_nesting
FUNCTION
Calculate the nesting function and output data to file.
INPUTS
ebands<ebands_t>=the ebands_t datatype cryst<crystal_t>=Info on unit cell and symmetries. filepath=File name for output data. prtnest = flags governing the format of the output file. see mknesting. 1 for X-Y format, 2 for XCrysden format (XSF) tsmear=Broadening used to approximation the delta function. fermie_nest qpath_vertices = vertices of the reciprocal space trajectory
OUTPUT
Return non-zero exist status if netsting factor cannot be produced. The errmsg string gives information on the error.
SIDE EFFECTS
Write data to file.
SOURCE
3704 integer function ebands_write_nesting(ebands,cryst,filepath,prtnest,tsmear,fermie_nest,qpath_vertices,errmsg) result(skip) 3705 3706 !Arguments ------------------------------------ 3707 class(ebands_t),intent(in) :: ebands 3708 type(crystal_t),intent(in) :: cryst 3709 integer,intent(in) :: prtnest 3710 real(dp),intent(in) :: tsmear,fermie_nest 3711 character(len=*),intent(in) :: filepath 3712 character(len=*),intent(out) :: errmsg 3713 !arrays 3714 real(dp),intent(in) :: qpath_vertices(:,:) 3715 3716 !Local variables------------------------------- 3717 !scalaras 3718 integer :: ikpt,spin,iband,nqpath 3719 real(dp) :: invgauwidth,prefact,fermie 3720 !arrays 3721 real(dp), allocatable :: fs_weights(:,:,:) 3722 3723 ! ********************************************************************* 3724 3725 skip = 0; errmsg = "" 3726 if (any(ebands%nband /= ebands%nband(1))) then 3727 errmsg = 'mknesting can not handle variable nband(1:nkpt). Skipped.'//& 3728 ch10//' Correct input file to get nesting output' 3729 skip = 1; return 3730 end if 3731 3732 if (ebands%nshiftk /= 1) then 3733 errmsg = 'mknesting does not support nshiftk > 1. Change ngkpt and shiftk to have only one shift after inkpts' 3734 skip = 1; return 3735 end if 3736 3737 ! FIXME: needs to be generalized to complete the k grid for one of the arguments to mknesting 3738 fermie = ebands%fermie 3739 nqpath = size(qpath_vertices, dim=2) 3740 3741 ! Compute weights. Set sigma to 0.1 eV is tsmear is zero 3742 invgauwidth = one / (0.1_dp * eV_Ha); if (tsmear > tol10) invgauwidth = one / tsmear 3743 prefact = one / sqrt(pi) * invgauwidth 3744 3745 ABI_MALLOC(fs_weights, (ebands%nband(1), ebands%nkpt, ebands%nsppol)) 3746 3747 do spin=1,ebands%nsppol 3748 do ikpt=1,ebands%nkpt 3749 do iband=1,ebands%nband(1) 3750 fs_weights(iband, ikpt, spin) = prefact * & 3751 exp(-(invgauwidth*(ebands%eig(iband,ikpt,spin)-(fermie + fermie_nest)))**2) 3752 end do 3753 end do 3754 end do 3755 3756 if (any(ebands%kptopt == [3, 4])) then ! no symmetry 3757 call mknesting(ebands%nkpt,ebands%kptns,ebands%kptrlatt,ebands%nband(1),fs_weights,nqpath,& 3758 qpath_vertices,1,[zero, zero, zero],filepath,cryst%gprimd,cryst%gmet,prtnest,identity_3d) 3759 else 3760 call mknesting(ebands%nkpt,ebands%kptns,ebands%kptrlatt,ebands%nband(1),fs_weights,nqpath,& 3761 qpath_vertices,1, [zero, zero, zero], filepath,cryst%gprimd,cryst%gmet,prtnest,identity_3d,& 3762 nsym=cryst%nsym, symrec=cryst%symrec) 3763 end if 3764 3765 ABI_FREE(fs_weights) 3766 3767 end function ebands_write_nesting
m_ebands/ebands_write_xmgrace [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
ebands_write_xmgrace
FUNCTION
Write bands in Xmgrace format. This routine should be called by a single processor. Use the driver `ebands_write` to support different formats.
INPUTS
filename=Filename [kptbounds(:,:)]=Optional argument giving the extrema of the k-path.
OUTPUT
Only writing
SOURCE
5570 subroutine ebands_write_xmgrace(ebands, filename, kptbounds) 5571 5572 !Arguments ------------------------------------ 5573 !scalars 5574 type(ebands_t),intent(in) :: ebands 5575 character(len=*),intent(in) :: filename 5576 real(dp),optional,intent(in) :: kptbounds(:,:) 5577 5578 !Local variables------------------------------- 5579 !scalars 5580 integer :: unt,ik,spin,band,ii,start,nkbounds 5581 character(len=500) :: msg 5582 !arrays 5583 integer :: g0(3) 5584 integer,allocatable :: bounds2kpt(:) 5585 5586 ! ********************************************************************* 5587 5588 nkbounds = 0 5589 if (present(kptbounds)) then 5590 if (product(shape(kptbounds)) > 0 ) then 5591 ! Find correspondence between kptbounds and k-points in ebands. 5592 nkbounds = size(kptbounds, dim=2) 5593 ABI_MALLOC(bounds2kpt, (nkbounds)) 5594 bounds2kpt = 1; start = 1 5595 do ii=1,nkbounds 5596 do ik=start,ebands%nkpt 5597 if (isamek(ebands%kptns(:, ik), kptbounds(:, ii), g0)) then 5598 bounds2kpt(ii) = ik; start = ik + 1; exit 5599 end if 5600 end do 5601 end do 5602 end if 5603 end if 5604 5605 if (open_file(filename, msg, newunit=unt, form="formatted", action="write") /= 0) then 5606 ABI_ERROR(msg) 5607 end if 5608 5609 write(unt,'(a)') "# Grace project file" 5610 write(unt,'(a)') "# Generated by Abinit" 5611 write(unt,'(4(a,i0))') & 5612 "# mband: ",ebands%mband,", nkpt: ",ebands%nkpt,", nsppol: ",ebands%nsppol,", nspinor: ",ebands%nspinor 5613 write(unt,'(a,f8.2,a,i0,2(a,f8.2))') & 5614 "# nelect: ",ebands%nelect,", occopt: ",ebands%occopt,", tsmear: ",ebands%tsmear,", tphysel: ",ebands%tphysel 5615 write(unt,'(a,f8.2,a)') "# Energies are in eV. Zero set to efermi, previously it was at: ",ebands%fermie * Ha_eV, " [eV]" 5616 write(unt,'(a)')"# List of k-points and their index (C notation i.e. count from 0)" 5617 do ik=1,ebands%nkpt 5618 write(unt, "(a)")sjoin("#", itoa(ik-1), ktoa(ebands%kptns(:,ik))) 5619 end do 5620 write(unt,'(a)') "@page size 792, 612" 5621 write(unt,'(a)') "@page scroll 5%" 5622 write(unt,'(a)') "@page inout 5%" 5623 write(unt,'(a)') "@link page off" 5624 write(unt,'(a)') "@with g0" 5625 write(unt,'(a)') "@world xmin 0.00" 5626 write(unt,'(a,i0)') '@world xmax ',ebands%nkpt 5627 write(unt,'(a,es16.8)') '@world ymin ',minval((ebands%eig - ebands%fermie) * Ha_eV) 5628 write(unt,'(a,es16.8)') '@world ymax ',maxval((ebands%eig - ebands%fermie) * Ha_eV) 5629 write(unt,'(a)') '@default linewidth 1.5' 5630 write(unt,'(a)') '@xaxis tick on' 5631 write(unt,'(a)') '@xaxis tick major 1' 5632 write(unt,'(a)') '@xaxis tick major color 1' 5633 write(unt,'(a)') '@xaxis tick major linestyle 3' 5634 write(unt,'(a)') '@xaxis tick major grid on' 5635 write(unt,'(a)') '@xaxis tick spec type both' 5636 write(unt,'(a)') '@xaxis tick major 0, 0' 5637 if (nkbounds /= 0) then 5638 write(unt,'(a,i0)') '@xaxis tick spec ',nkbounds 5639 do ik=1,nkbounds 5640 !write(unt,'(a,i0,a,a)') '@xaxis ticklabel ',ik-1,',', "foo" 5641 write(unt,'(a,i0,a,i0)') '@xaxis tick major ',ik-1,' , ',bounds2kpt(ik) - 1 5642 end do 5643 end if 5644 write(unt,'(a)') '@xaxis ticklabel char size 1.500000' 5645 write(unt,'(a)') '@yaxis tick major 10' 5646 write(unt,'(a)') '@yaxis label "Band Energy [eV]"' 5647 write(unt,'(a)') '@yaxis label char size 1.500000' 5648 write(unt,'(a)') '@yaxis ticklabel char size 1.500000' 5649 ii = -1 5650 do spin=1,ebands%nsppol 5651 do band=1,ebands%mband 5652 ii = ii + 1 5653 write(unt,'(a,i0,a,i0)') '@ s',ii,' line color ',spin 5654 end do 5655 end do 5656 ii = -1 5657 do spin=1,ebands%nsppol 5658 do band=1,ebands%mband 5659 ii = ii + 1 5660 write(unt,'(a,i0)') '@target G0.S',ii 5661 write(unt,'(a)') '@type xy' 5662 do ik=1,ebands%nkpt 5663 write(unt,'(i0,1x,es16.8)') ik-1, (ebands%eig(band, ik, spin) - ebands%fermie) * Ha_eV 5664 end do 5665 write(unt,'(a)') '&' 5666 end do 5667 end do 5668 5669 close(unt) 5670 5671 ABI_SFREE(bounds2kpt) 5672 5673 end subroutine ebands_write_xmgrace
m_ebands/edos_free [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
edos_free
FUNCTION
Free the memory allocated in edos_t
SOURCE
3331 subroutine edos_free(edos) 3332 3333 !Arguments ------------------------------------ 3334 class(edos_t),intent(inout) :: edos 3335 3336 ! ********************************************************************* 3337 3338 !real 3339 ABI_SFREE(edos%mesh) 3340 ABI_SFREE(edos%dos) 3341 ABI_SFREE(edos%idos) 3342 ABI_SFREE(edos%gef) 3343 ABI_SFREE(edos%ghf) 3344 3345 end subroutine edos_free
m_ebands/edos_get_carriers [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
edos_get_carriers
FUNCTION
Compute number of holes (nh) and electrons (ne) per unit cell from a given list of `ntemp` temperatures `kTmesh` and chemical potentials `mu_e`. Return n_ehst(2, nsppol, ntemp) where the first dimension if for electrons/holes. If nsppol == 2, the second dimension is the number of e/h for spin else the total number of e/h summed over spins.!! To discern between electrons and holes in semiconductors we assume that ef is inside the gap.
INPUTS
OUTPUT
SOURCE
3623 subroutine edos_get_carriers(edos, ntemp, kTmesh, mu_e, n_ehst) 3624 3625 !Arguments ------------------------------------ 3626 class(edos_t),intent(in) :: edos 3627 integer,intent(in) :: ntemp 3628 !arrays 3629 real(dp),intent(in) :: kTmesh(ntemp), mu_e(ntemp) 3630 real(dp),intent(out) :: n_ehst(2, edos%nsppol, ntemp) 3631 3632 !Local variables------------------------------- 3633 integer :: itemp, iw, spin 3634 real(dp),allocatable :: values(:) 3635 3636 ! ************************************************************************* 3637 3638 ! Copy important dimensions 3639 n_ehst = zero 3640 ABI_MALLOC(values, (edos%nw)) 3641 3642 do itemp=1,ntemp 3643 3644 ! For electrons (assuming ef inside the gap if semiconductor) 3645 do spin=1,edos%nsppol 3646 do iw=1,edos%nw 3647 if (edos%mesh(iw) >= mu_e(itemp)) then 3648 values(iw) = edos%dos(iw, spin) * occ_fd(edos%mesh(iw), kTmesh(itemp), mu_e(itemp)) 3649 else 3650 values(iw) = zero 3651 end if 3652 end do 3653 n_ehst(1, spin, itemp) = simpson(edos%step, values) 3654 end do ! spin 3655 3656 ! For holes 3657 do spin=1,edos%nsppol 3658 do iw=1,edos%nw 3659 if (edos%mesh(iw) < mu_e(itemp)) then 3660 values(iw) = edos%dos(iw, spin) * (one - occ_fd(edos%mesh(iw), kTmesh(itemp), mu_e(itemp))) 3661 else 3662 values(iw) = zero 3663 end if 3664 end do 3665 n_ehst(2, spin, itemp) = simpson(edos%step, values) 3666 end do ! spin 3667 3668 end do ! itemp 3669 3670 if (edos%nsppol == 1 .and. edos%nspinor == 1) n_ehst = two * n_ehst 3671 ABI_FREE(values) 3672 3673 end subroutine edos_get_carriers
m_ebands/edos_ncwrite [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
edos_ncwrite
FUNCTION
Write results to netcdf file.
INPUTS
edos<edos_t>=DOS container ncid=NC file handle. [prefix]=String prepended to netcdf dimensions/variables (HDF5 poor-man groups) Empty string if not specified.
OUTPUT
ncerr= netcdf exit status.
SOURCE
3463 integer function edos_ncwrite(edos, ncid, prefix) result(ncerr) 3464 3465 !Arguments ------------------------------------ 3466 integer,intent(in) :: ncid 3467 class(edos_t),intent(in) :: edos 3468 character(len=*),optional,intent(in) :: prefix 3469 3470 !Local variables------------------------------- 3471 character(len=500) :: prefix_ 3472 3473 ! ************************************************************************* 3474 3475 prefix_ = ""; if (present(prefix)) prefix_ = trim(prefix) 3476 3477 ! Define dimensions. 3478 ncerr = nctk_def_dims(ncid, [ & 3479 nctkdim_t("nsppol_plus1", edos%nsppol + 1), nctkdim_t("edos_nw", edos%nw)], defmode=.True., prefix=prefix_) 3480 NCF_CHECK(ncerr) 3481 3482 ! Define variables 3483 NCF_CHECK(nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "edos_intmeth", "edos_nkibz"], prefix=prefix_)) 3484 NCF_CHECK(nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "edos_ief", "edos_ihf"], prefix=prefix_)) 3485 NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "edos_broad"], prefix=prefix_)) 3486 3487 ncerr = nctk_def_arrays(ncid, [ & 3488 nctkarr_t("edos_mesh", "dp", "edos_nw"), & 3489 nctkarr_t("edos_dos", "dp", "edos_nw, nsppol_plus1"), & 3490 nctkarr_t("edos_idos", "dp", "edos_nw, nsppol_plus1"), & 3491 nctkarr_t("edos_gef", "dp", "nsppol_plus1"), & 3492 nctkarr_t("edos_ghf", "dp", "nsppol_plus1") & 3493 ], prefix=prefix_) 3494 NCF_CHECK(ncerr) 3495 3496 ! Write data. 3497 NCF_CHECK(nctk_set_datamode(ncid)) 3498 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_intmeth")), edos%intmeth)) 3499 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_nkibz")), edos%nkibz)) 3500 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_ief")), edos%ief)) 3501 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_ihf")), edos%ihf)) 3502 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_broad")), edos%broad)) 3503 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_mesh")), edos%mesh)) 3504 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_dos")), edos%dos)) 3505 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_idos")), edos%idos)) 3506 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_gef")), edos%gef)) 3507 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_ghf")), edos%ghf)) 3508 3509 contains 3510 pure function pre(istr) result(ostr) 3511 character(len=*),intent(in) :: istr 3512 character(len=len_trim(prefix_) + len_trim(istr)+1) :: ostr 3513 ostr = trim(prefix_) // trim(istr) 3514 end function pre 3515 3516 end function edos_ncwrite
m_ebands/edos_print [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
edos_print
FUNCTION
Print DOS info to Fortran unit.
INPUTS
edos<edos_t>=DOS container [unit]=Unit number for output. Defaults to std_out
OUTPUT
Only writing.
SOURCE
3535 subroutine edos_print(edos, unit, header) 3536 3537 !Arguments ------------------------------------ 3538 class(edos_t),intent(in) :: edos 3539 integer,optional,intent(in) :: unit 3540 character(len=*),optional,intent(in) :: header 3541 3542 !Local variables------------------------------- 3543 integer :: unt 3544 3545 ! ************************************************************************* 3546 3547 unt = std_out; if (present(unit)) unt = unit 3548 if (unt == dev_null) return 3549 3550 if (present(header)) then 3551 write(unt, "(a)") ch10//' === '//trim(adjustl(header))//' === ' 3552 else 3553 write(unt, "(a)") ch10 3554 end if 3555 3556 select case (edos%intmeth) 3557 case (1) 3558 write(unt, "(a,f5.1,a)") " Gaussian method with broadening: ", edos%broad * Ha_meV, " (meV)" 3559 case (2) 3560 write(unt, "(a)")" Linear tetrahedron method." 3561 case (-2) 3562 write(unt, "(a)")" Linear tetrahedron method with Blochl corrections." 3563 case default 3564 ABI_ERROR(sjoin("Wrong intmeth:", itoa(edos%intmeth))) 3565 end select 3566 3567 write(unt, "(a,f5.1,a, i0)")" Mesh step: ", edos%step * Ha_meV, " (meV) with npts: ", edos%nw 3568 write(unt, "(2(a,f5.1),a)")" From emin: ", edos%mesh(1) * Ha_eV, " to emax: ", edos%mesh(edos%nw) * Ha_eV, " (eV)" 3569 write(unt, "(a, i0)")" Number of k-points in the IBZ: ", edos%nkibz 3570 3571 if (edos%ief == 0) then 3572 write(unt, "(a, /)")" edos%ief == 0 --> Cannot print quantities at the Fermi level." 3573 return 3574 end if 3575 3576 write(unt,'(a,es16.8,a)')' Fermi level: ',edos%mesh(edos%ief) * Ha_eV, " (eV)" 3577 write(unt,"(a,es16.8)")" Total electron DOS at Fermi level in states/eV: ", edos%gef(0) / Ha_eV 3578 if (edos%nsppol == 2) then 3579 write(unt,"(a,es16.8)")" g(eF) for spin up: ", edos%gef(1) / Ha_eV 3580 write(unt,"(a,es16.8)")" g(eF) for spin down:", edos%gef(2) / Ha_eV 3581 end if 3582 write(unt,"(a,f6.1)")" Total number of electrons at eF: ", edos%idos(edos%ief, 0) 3583 if (edos%nsppol == 2) then 3584 write(unt,"(a,es16.8)")" IDOS(eF) for spin up: ", edos%idos(edos%ief, 1) 3585 write(unt,"(a,es16.8)")" IDOS(eF) for spin down:", edos%idos(edos%ief, 2) 3586 end if 3587 if (edos%ihf /= edos%ief) then 3588 write(unt,'(a,es16.8,a)')' Fermi level for excited holes: ',edos%mesh(edos%ihf) * Ha_eV, " (eV)" 3589 write(unt,"(a,es16.8)")" Total hole DOS at Fermi level in states/eV: ", edos%ghf(0) / Ha_eV 3590 if (edos%nsppol == 2) then 3591 write(unt,"(a,es16.8)")" g(hF) for spin up: ", edos%ghf(1) / Ha_eV 3592 write(unt,"(a,es16.8)")" g(hF) for spin down:", edos%ghf(2) / Ha_eV 3593 end if 3594 write(unt,"(a,f6.1)")" Total number of electrons at hF: ", edos%idos(edos%ihf, 0) 3595 if (edos%nsppol == 2) then 3596 write(unt,"(a,es16.8)")" N(hF) for spin up: ", edos%idos(edos%ihf, 1) 3597 write(unt,"(a,es16.8)")" N(hF) for spin down:", edos%idos(edos%ihf, 2) 3598 end if 3599 end if 3600 3601 write(unt, "(a)")"" 3602 3603 end subroutine edos_print
m_ebands/edos_t [ Types ]
[ Top ] [ m_ebands ] [ Types ]
NAME
edos_t
FUNCTION
Store the electronic DOS
SOURCE
133 type,public :: edos_t 134 135 integer :: nsppol = -1 136 ! Number of spins. 137 138 integer :: nspinor = -1 139 ! Number of spinors 140 141 integer :: nkibz = -1 142 ! Number of k-points in the IBZ. 143 144 integer :: nw = -1 145 ! Number of points in the frequency mesh. 146 147 integer :: ief = 0 148 ! Rightmost Index of the energy mesh such as IDOS[mesh[ief]] < nelect. 149 ! 0 if Fermi level could not be computed 150 ! Note the value of gef stored in edos_t is computed by performing 151 ! a linear interpolation between ief and ief + 1 152 153 integer :: ihf = 0 154 ! Like ief (see above) 155 ! But for fermi level of thermalized holes in valence bands 156 ! For occopt 9 purposes 157 158 integer :: intmeth = 0 159 ! 1 for gaussian, 2 tetra 160 161 real(dp) :: broad = zero 162 ! Gaussian broadening 163 164 real(dp) :: step = -one 165 ! Step of the mesh 166 167 real(dp) :: nelect = zero 168 ! Number of electrons taken from ebands. 169 170 real(dp),allocatable :: mesh(:) 171 ! mesh(nw) 172 173 real(dp),allocatable :: dos(:,:) 174 ! dos(nw, 0:nsppol) 175 ! Total DOS, spin up and spin down component. 176 177 real(dp),allocatable :: idos(:,:) 178 ! idos(nw, 0:nsppol) 179 ! Integrated DOS: (total, spin up, spin down) component. 180 181 real(dp),allocatable :: gef(:) 182 ! gef(0:nsppol) 183 ! DOS at the Fermi level. Total, spin up, spin down 184 185 real(dp),allocatable :: ghf(:) 186 ! ghf(0:nsppol) 187 ! DOS at the Fermi level of thermalized holes. Total, spin up, spin down 188 189 contains 190 191 procedure :: free => edos_free 192 ! Free memory 193 194 procedure :: write => edos_write 195 ! Write results to file (formatted mode) 196 197 procedure :: print => edos_print 198 ! Print eDOS info to Fortran unit. 199 200 procedure :: ncwrite => edos_ncwrite 201 ! Write eDOS to netcdf file. 202 203 procedure :: get_carriers => edos_get_carriers 204 ! Compute number of holes (nh) and electrons (ne) per unit cell from a given 205 ! list of `ntemp` temperatures `kTmesh` and chemical potentials `mu_e`. 206 207 end type edos_t
m_ebands/edos_write [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
edos_write
FUNCTION
Write results to file.
INPUTS
edos<edos_t>=DOS container path=File name.
OUTPUT
Only writing.
SOURCE
3366 subroutine edos_write(edos, path) 3367 3368 !Arguments ------------------------------------ 3369 character(len=*),intent(in) :: path 3370 class(edos_t),intent(in) :: edos 3371 3372 !Local variables------------------------------- 3373 integer :: iw,spin,unt 3374 real(dp) :: efermi, gef_tot, gef_up, gef_down 3375 character(len=500) :: msg 3376 type(yamldoc_t) :: ydoc 3377 3378 ! ************************************************************************* 3379 3380 if (open_file(path, msg, newunit=unt, form="formatted", action="write") /= 0) then 3381 ABI_ERROR(msg) 3382 end if 3383 3384 ! Write header (human-readable format) 3385 write(unt,'(a)')'# Electron density of states: Energy in eV, DOS in states/eV per unit cell.' 3386 write(unt,"(a)")"# The zero of energies corresponds to the Fermi level." 3387 3388 ! Add Yaml section with parameters. 3389 ydoc = yamldoc_open("EDOS_PARAMS") 3390 call ydoc%add_ints("nkibz, nsppol, nspinor, intmeth, edos_npts", & 3391 [edos%nkibz, edos%nsppol, edos%nspinor, edos%intmeth, edos%nw]) 3392 call ydoc%add_reals("nelect, edos_mesh_step_eV", & 3393 [edos%nelect, edos%step * Ha_eV]) 3394 3395 select case (edos%intmeth) 3396 case (1) 3397 call ydoc%add_string("method", "Gaussian") 3398 call ydoc%add_real("gaussian_broadening_eV", edos%broad * Ha_eV) 3399 case (2) 3400 call ydoc%add_string("method", "Linear tetrahedron") 3401 case (-2) 3402 call ydoc%add_string("method", "Linear tetrahedron method with Blochl corrections") 3403 case default 3404 ABI_ERROR(sjoin("Wrong method:", itoa(edos%intmeth))) 3405 end select 3406 3407 if (edos%ief == 0) then 3408 call ydoc%set_keys_to_string("Fermi_level_eV, gef, gef_up, gef_down", "null") 3409 efermi = zero 3410 else 3411 efermi = edos%mesh(edos%ief) 3412 call ydoc%add_real("Fermi_level_eV", efermi * Ha_eV) 3413 gef_tot = edos%gef(0) / Ha_eV 3414 gef_up = gef_tot / two; gef_down = gef_tot / two 3415 if (edos%nsppol == 2) then 3416 gef_up = edos%gef(1) / Ha_eV; gef_down = edos%gef(2) / Ha_eV 3417 end if 3418 if (edos%nspinor == 1) then 3419 call ydoc%add_reals("gef, gef_up, gef_down", [gef_tot, gef_up, gef_down]) 3420 else 3421 call ydoc%add_reals("gef", [gef_tot]) 3422 end if 3423 end if 3424 3425 ! Write header in Yaml format but prepend # so that one can still use tools such as gnuplot or xmgrace. 3426 call ydoc%write_and_free(unt, firstchar="#") 3427 3428 ! Write data. 3429 write(unt,"(a)")"# Energy DOS_TOT IDOS_TOT DOS[spin=UP] IDOS[spin=UP] ..." 3430 do iw=1,edos%nw 3431 write(unt,'(es17.8)',advance='no')(edos%mesh(iw) - efermi) * Ha_eV 3432 do spin=0,edos%nsppol 3433 write(unt,'(2es17.8)',advance='no')max(edos%dos(iw,spin) / Ha_eV, tol30), max(edos%idos(iw,spin), tol30) 3434 end do 3435 write(unt,*) 3436 end do 3437 3438 close(unt) 3439 3440 end subroutine edos_write
m_ebands/gaps_free [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
gaps_free
FUNCTION
Free the memory allocated in gaps_t
SOURCE
617 subroutine gaps_free(gaps) 618 619 !Arguments ------------------------------------ 620 class(gaps_t),intent(inout) :: gaps 621 622 ! ********************************************************************* 623 624 !integer 625 ABI_SFREE(gaps%fo_kpos) 626 ABI_SFREE(gaps%ierr) 627 628 !real 629 ABI_SFREE(gaps%fo_values) 630 ABI_SFREE(gaps%vb_max) 631 ABI_SFREE(gaps%cb_min) 632 ABI_SFREE(gaps%optical_kpoints) 633 ABI_SFREE(gaps%fund_kpoints) 634 635 !chars 636 ABI_SFREE(gaps%errmsg_spin) 637 638 end subroutine gaps_free
m_ebands/gaps_print [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
gaps_print
FUNCTION
Print info on the fundamental and direct gap.
INPUTS
gaps<gaps_t>=Object with info on the gaps. [header]=Optional title. [unit]=Optional unit for output (std_out if not specified) [kTmesh]=List of temperatures. If present activates output of (T, mu_e, band_edges) [mu_e]=List of Fermi levels for each T.
OUTPUT
Only writing.
SOURCE
662 subroutine gaps_print(gaps, unit, header, kTmesh, mu_e) 663 664 !Arguments ------------------------------------ 665 !scalars 666 class(gaps_t),intent(in) :: gaps 667 integer,intent(in),optional :: unit 668 character(len=*),intent(in),optional :: header 669 real(dp),optional,intent(in) :: kTmesh(:), mu_e(:) 670 671 !Local variables------------------------------- 672 !scalars 673 integer :: spin, ikopt, ivk, ick, unt, itemp, ntemp 674 real(dp) :: fun_gap, opt_gap, csi_c, csi_v 675 character(len=500) :: msg 676 677 ! ********************************************************************* 678 679 unt = std_out; if (present(unit)) unt = unit 680 if (unt == dev_null) return 681 682 do spin=1,gaps%nsppol 683 if (spin == 1) then 684 msg = ch10 685 if (present(header)) msg = ch10//' === '//trim(adjustl(header))//' === ' 686 call wrtout(unt, msg) 687 end if 688 689 if (gaps%ierr(spin) /= 0) then 690 call wrtout(unt, gaps%errmsg_spin(spin)) 691 continue 692 end if 693 694 ! Get minimum of the direct Gap. 695 fun_gap = gaps%fo_values(1, spin) 696 opt_gap = gaps%fo_values(2, spin) 697 698 if (any(gaps%fo_kpos(:,spin) == 0)) then 699 call wrtout(unt, sjoin(" Cannot detect gap for spin: ", itoa(spin))) 700 cycle 701 end if 702 703 ivk = gaps%fo_kpos(1, spin) 704 ick = gaps%fo_kpos(2, spin) 705 ikopt = gaps%fo_kpos(3, spin) 706 707 ! >>>> For spin 2 708 !Direct band gap semiconductor. 709 !Fundamental gap: 4.48 (eV) 710 ! VBM: 4.47 (eV) at k: [ 0.0000E+00, 0.0000E+00, 0.0000E+00] 711 ! CBM: 8.96 (eV) at k: [ 0.0000E+00, 0.0000E+00, 0.0000E+00] 712 !Optical gap: 4.48 (eV) at k:[ 0.0000E+00, 0.0000E+00, 0.0000E+00] 713 714 if (gaps%nsppol == 2) call wrtout(unt, sjoin(' >>>> For spin ', itoa(spin))) 715 if (ivk == ick) call wrtout(unt, " Direct band gap semiconductor") 716 if (ivk /= ick) call wrtout(unt, " Indirect band gap semiconductor") 717 write(msg, "(a,f9.3,a )")" Fundamental gap: ", fun_gap * Ha_eV, " (eV)" 718 call wrtout(unt, msg) 719 write(msg, "(a,f9.3,2a)")" VBM: ", gaps%vb_max(spin) * Ha_eV, " (eV) at k: ", trim(ktoa(gaps%fund_kpoints(:,1,spin))) 720 call wrtout(unt, msg) 721 write(msg, "(a,f9.3,2a)")" CBM: ", gaps%cb_min(spin) * Ha_eV, " (eV) at k: ", trim(ktoa(gaps%fund_kpoints(:,2,spin))) 722 call wrtout(unt, msg) 723 write(msg, "(a,f9.3,2a)")" Direct gap: ", opt_gap * Ha_eV," (eV) at k: ", trim(ktoa(gaps%optical_kpoints(:,spin))) 724 call wrtout(unt, msg) 725 !write(msg, "((2(a, f9.3)))")" Fermi level:", gaps%fermie * Ha_eV, " (eV) with nelect:", gaps%nelect 726 !call wrtout(unt, msg) 727 728 if (present(mu_e) .and. present(kTmesh) .and. all(gaps%ierr == 0)) then 729 ntemp = size(mu_e) 730 call wrtout(unt, " Position of CBM/VBM with respect to the Fermi level:", pre_newlines=1) 731 call wrtout(unt, " Notations: mu_e = Fermi level, D_v = (mu_e - VBM), D_c = (CBM - mu_e)") 732 call wrtout(unt, " T(K) kT (eV) mu_e (eV) D_v (eV) D_c (eV)", pre_newlines=1) 733 do itemp=1,ntemp 734 csi_c = gaps%cb_min(spin) - mu_e(itemp) 735 csi_v = -gaps%vb_max(spin) + mu_e(itemp) 736 write(msg, "(f6.1, 1x, 4(f9.3, 1x))") & 737 kTmesh(itemp) / kb_HaK, kTmesh(itemp) * Ha_eV, mu_e(itemp) * Ha_eV, csi_v * Ha_eV, csi_c * Ha_eV 738 call wrtout(unt, msg) 739 end do 740 call wrtout(unt, "") 741 end if 742 743 end do ! spin 744 745 if (any(gaps%fo_kpos == 0)) then 746 write(msg, "((2(a, f9.3)))") " Fermi level:", gaps%fermie * Ha_eV, " (eV) with nelect:", gaps%nelect 747 call wrtout(unt, msg) 748 end if 749 750 call wrtout(unt, "") 751 752 end subroutine gaps_print
m_ebands/gaps_t [ Types ]
[ Top ] [ m_ebands ] [ Types ]
NAME
gaps_t
FUNCTION
Structure with information on the fundamental and direct gaps returned by ebands_report_gap.
SOURCE
274 type,public :: gaps_t 275 276 integer :: nsppol 277 ! Number of spins. 278 279 integer,allocatable :: fo_kpos(:,:) 280 ! fo_kpos(3,nsppol) 281 ! fo_kpos(1:2,spin) ==> Indices of the k-points where the homo, lumo states are located (for each spin). 282 ! fo_kpos(3,spin) ==> the index of k-point where the direct gap is located (for each spin). 283 284 real(dp) :: fermie 285 ! Fermi energy taken from ebands. 286 287 real(dp) :: nelect 288 ! Number of electrons taken from ebands. 289 290 integer,allocatable :: ierr(:) 291 ! ierr(nsppol 292 ! 0 if the gap has been computed. 293 ! 1 if the system (or spin-channel) is metallic. 294 ! 2 if gaps were not computed (because there are only valence bands). 295 296 real(dp),allocatable :: fo_values(:,:) 297 ! fo_values(2,nsppol)] 298 ! Fundamental and direct gaps (in Hartree) for each spin. 299 300 real(dp),allocatable :: vb_max(:), cb_min(:) 301 ! vb_max(nsppol) 302 ! valence band max and conduction band min for each spin in Ha. 303 ! Only for Semiconductors, set to (+, -) huge(one) for metals. 304 305 real(dp),allocatable :: optical_kpoints(:,:) 306 ! (3, nsppol) 307 ! kpoint of optical gap for each spin 308 309 real(dp),allocatable :: fund_kpoints(:,:, :) 310 ! (3, 2, nsppol) 311 ! kpoint of the fundamental gap for (val, cond) and each spin 312 313 character(len=500),allocatable :: errmsg_spin(:) 314 ! errmsg_spin(nsppol) 315 ! String with human-readable error messages if ierr(spin) != 0. 316 317 contains 318 319 procedure :: free => gaps_free 320 ! Free memory 321 322 procedure :: print => gaps_print 323 ! Print info on the gaps 324 325 end type gaps_t 326 327 public :: ebands_get_gaps ! Build the gaps object from a bandstructure. 328 public :: ebands_print_gaps ! Helper function to print gaps directrly from ebands.
m_ebands/get_eneocc_vect [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
get_eneocc_vect
FUNCTION
Retrieve energies or occupations from a ebands_t structure accessing by name. Results are reported in a vector to facilitate the interface with other abinit routines.
INPUTS
ebands<ebands_t>The type containing the data. arr_name=The name of the quantity to retrieve. Allowed values are == "eig" == For the eigenvalues. == "occ" == For the occupation numbers. == "doccde" == For the derivative of the occupancies wrt the energy.
OUTPUT
vect(ebands%bantot)=The values required.
SOURCE
1374 subroutine get_eneocc_vect(ebands, arr_name, vect) 1375 1376 !Arguments ------------------------------------ 1377 !scalars 1378 character(len=*),intent(in) :: arr_name 1379 class(ebands_t),intent(in) :: ebands 1380 real(dp),intent(out) :: vect(ebands%bantot) 1381 1382 !Local variables------------------------------- 1383 integer :: nkpt,nsppol,mband,bantot 1384 ! ************************************************************************* 1385 1386 mband = ebands%mband; bantot = ebands%bantot; nkpt = ebands%nkpt; nsppol = ebands%nsppol 1387 1388 select case (arr_name) 1389 case ('occ') 1390 call pack_eneocc(nkpt, nsppol, mband, ebands%nband, bantot, ebands%occ, vect) 1391 case ('eig') 1392 call pack_eneocc(nkpt,nsppol,mband,ebands%nband,bantot,ebands%eig, vect) 1393 case ('doccde') 1394 call pack_eneocc(nkpt,nsppol,mband,ebands%nband,bantot,ebands%doccde,vect) 1395 case default 1396 ABI_BUG(sjoin('Wrong arr_name:', arr_name)) 1397 end select 1398 1399 end subroutine get_eneocc_vect
m_ebands/get_gaps_ [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
get_gaps_
FUNCTION
Private function that returns a gaps_t object with info on the fundamental and direct gap.
INPUTS
ebands<ebands_t>=Info on the band structure, the smearing technique and the physical temperature used.
OUTPUT
ierr=Return code (!=0 signals failure) gaps<gaps_t>=object with info on the gaps (caller is responsible for freeing the object).
SOURCE
500 type(gaps_t) function get_gaps_(ebands, ierr) result(gaps) 501 502 !Arguments ------------------------------------ 503 !scalars 504 class(ebands_t),target,intent(in) :: ebands 505 integer,intent(out) :: ierr 506 507 !Local variables------------------------------- 508 !scalars 509 integer :: ikibz,nband_k,spin,nsppol,ikopt,ivk,ick,ivb,icb 510 real(dp),parameter :: tol_fermi = tol6 511 real(dp) :: fun_gap, opt_gap 512 logical :: ismetal 513 !type(ebands_t) :: tmp_ebands 514 !arrays 515 integer :: val_idx(ebands%nkpt, ebands%nsppol) 516 real(dp) :: top_valence(ebands%nkpt), bot_conduct(ebands%nkpt) 517 518 ! ********************************************************************* 519 520 nsppol = ebands%nsppol 521 522 ! Initialize gaps_t 523 gaps%nsppol = nsppol 524 gaps%nelect = ebands%nelect 525 ABI_MALLOC(gaps%fo_kpos, (3, nsppol)) 526 ABI_MALLOC(gaps%ierr, (nsppol)) 527 ABI_MALLOC(gaps%fo_values, (2, nsppol)) 528 ABI_MALLOC(gaps%vb_max, (nsppol)) 529 ABI_MALLOC(gaps%cb_min, (nsppol)) 530 ABI_MALLOC(gaps%errmsg_spin, (nsppol)) 531 532 ABI_MALLOC(gaps%fund_kpoints, (3, 2, nsppol)) 533 ABI_MALLOC(gaps%optical_kpoints, (3, nsppol)) 534 gaps%fund_kpoints = huge(one) 535 gaps%optical_kpoints = huge(one) 536 537 gaps%fo_kpos = 0 538 gaps%ierr = 0 539 gaps%fo_values = zero 540 gaps%vb_max = huge(one); gaps%cb_min = -huge(one) 541 gaps%errmsg_spin(:) = "" 542 gaps%fermie = ebands%fermie 543 544 ! Compute "valence index" using efermi 545 val_idx(:,:) = ebands_get_valence_idx(ebands, tol_fermi=tol_fermi) 546 547 spin_loop: & 548 & do spin=1,nsppol 549 550 ! No output if system is metallic 551 ismetal = ANY(val_idx(:,spin) /= val_idx(1,spin)) 552 if (ismetal) then 553 gaps%ierr(spin) = 1 554 write(gaps%errmsg_spin(spin), "(a,i0)")" Detected metallic system for spin channel: ", spin 555 cycle 556 endif 557 558 ivb = val_idx(1, spin) 559 icb = ivb + 1 560 561 do ikibz=1,ebands%nkpt 562 nband_k = ebands%nband(ikibz + (spin-1)*ebands%nkpt) 563 top_valence(ikibz) = ebands%eig(ivb, ikibz, spin) 564 if (icb > nband_k) then 565 gaps%ierr(spin) = 2 566 gaps%errmsg_spin(spin) = "Not enough states to calculate the band gap." 567 cycle spin_loop 568 end if 569 bot_conduct(ikibz) = ebands%eig(icb, ikibz, spin) 570 end do 571 572 ! Minimum of the direct Gaps 573 ikopt = imin_loc(bot_conduct - top_valence) 574 opt_gap = bot_conduct(ikopt) - top_valence(ikopt) 575 576 ! Fundamental Gap 577 ick = imin_loc(bot_conduct) 578 ivk = imax_loc(top_valence) 579 580 gaps%vb_max(spin) = ebands%eig(ivb, ivk, spin) 581 gaps%cb_min(spin) = ebands%eig(icb, ick, spin) 582 fun_gap = ebands%eig(icb, ick, spin) - ebands%eig(ivb, ivk, spin) 583 gaps%fo_values(:, spin) = [fun_gap, opt_gap] 584 gaps%fo_kpos(:, spin) = [ivk, ick, ikopt] 585 586 gaps%optical_kpoints(:, spin) = ebands%kptns(:, ikopt) 587 gaps%fund_kpoints(:, 1, spin) = ebands%kptns(:, ivk) 588 gaps%fund_kpoints(:, 2, spin) = ebands%kptns(:, ick) 589 end do spin_loop 590 591 ierr = maxval(gaps%ierr) 592 593 if (ierr /= 0) then 594 ! Set VBM and CBM to fermie if metal. 595 do spin=1,nsppol 596 if (gaps%ierr(spin) /= 0) then 597 gaps%vb_max(spin) = ebands%fermie 598 gaps%cb_min(spin) = ebands%fermie 599 end if 600 end do 601 end if 602 603 end function get_gaps_
m_ebands/jdos_free [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
jdos_free
FUNCTION
Free memory
SOURCE
5078 subroutine jdos_free(jdos) 5079 5080 !Arguments ------------------------------------ 5081 class(jdos_t),intent(inout) :: jdos 5082 5083 ! ********************************************************************* 5084 5085 ABI_SFREE(jdos%mesh) 5086 ABI_SFREE(jdos%values) 5087 5088 end subroutine jdos_free
m_ebands/jdos_ncwrite [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
jdos_ncwrite
FUNCTION
Write JDOS to netcdf file.
INPUTS
ncid=NC file handle. [prefix]=String prepended to netcdf dimensions/variables (HDF5 poor-man groups) Empty string if not specified.
SOURCE
5019 integer function jdos_ncwrite(jdos, ncid, prefix) result(ncerr) 5020 5021 !Arguments ------------------------------------ 5022 !scalars 5023 class(jdos_t),intent(inout) :: jdos 5024 integer,intent(in) :: ncid 5025 character(len=*),optional,intent(in) :: prefix 5026 5027 !Local variables------------------------------- 5028 character(len=500) :: prefix_ 5029 5030 ! ********************************************************************* 5031 5032 prefix_ = ""; if (present(prefix)) prefix_ = trim(prefix) 5033 5034 ! Define dimensions. 5035 ncerr = nctk_def_dims(ncid, [ & 5036 nctkdim_t("nsppol_plus1", jdos%nsppol + 1), nctkdim_t("jdos_nw", jdos%nw)], defmode=.True., prefix=prefix_) 5037 NCF_CHECK(ncerr) 5038 5039 ! Define variables 5040 NCF_CHECK(nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "jdos_intmeth", "jdos_nkibz"], prefix=prefix_)) 5041 NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "jdos_broad"], prefix=prefix_)) 5042 5043 ncerr = nctk_def_arrays(ncid, [ & 5044 nctkarr_t("jdos_mesh", "dp", "jdos_nw"), & 5045 nctkarr_t("jdos_values", "dp", "jdos_nw, nsppol_plus1") & 5046 ], prefix=prefix_) 5047 NCF_CHECK(ncerr) 5048 5049 ! Write data. 5050 NCF_CHECK(nctk_set_datamode(ncid)) 5051 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("jdos_intmeth")), jdos%intmeth)) 5052 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("jdos_nkibz")), jdos%nkibz)) 5053 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("jdos_broad")), jdos%broad)) 5054 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("jdos_mesh")), jdos%mesh)) 5055 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("edos_values")), jdos%values)) 5056 5057 contains 5058 pure function pre(istr) result(ostr) 5059 character(len=*),intent(in) :: istr 5060 character(len=len_trim(prefix_) + len_trim(istr)+1) :: ostr 5061 ostr = trim(prefix_) // trim(istr) 5062 end function pre 5063 5064 end function jdos_ncwrite
m_ebands/jdos_t [ Types ]
[ Top ] [ m_ebands ] [ Types ]
NAME
jdos_t
FUNCTION
Store the electron joint DOS
SOURCE
219 type,public :: jdos_t 220 221 integer :: nsppol 222 ! Number of spins. 223 224 integer :: nkibz 225 ! Number of k-points in the IBZ. 226 227 integer :: nw 228 ! Number of points in the frequency mesh. 229 230 integer :: intmeth 231 ! 1 for gaussian, 2 tetra 232 233 real(dp) :: broad = zero 234 ! Gaussian broadening 235 236 real(dp) :: step 237 ! Step of the mesh 238 239 real(dp),allocatable :: mesh(:) 240 ! mesh(nw) 241 242 real(dp),allocatable :: values(:,:) 243 ! dos(nw,0:nsppol) 244 ! Total jDOS, spin up and spin down component. 245 246 contains 247 248 procedure :: free => jdos_free 249 ! Free memory 250 251 !procedure :: write => jdos_write 252 ! Write results to file (formatted mode) 253 254 !procedure :: print => jdos_print 255 ! Print jDOS info to Fortran unit. 256 257 procedure :: ncwrite => jdos_ncwrite 258 ! Write jDOS to netcdf file. 259 260 end type jdos_t
m_ebands/klinterp_eval_bsd [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
klinterp_eval_bsd
FUNCTION
INPUTS
OUTPUT
SOURCE
6047 subroutine klinterp_eval_bsd(self, kpt, vals_bsd) 6048 6049 !Arguments ------------------------------------ 6050 !scalars 6051 class(klinterp_t),intent(in) :: self 6052 real(dp),intent(in) :: kpt(3) 6053 real(dp),intent(out) :: vals_bsd(self%bsize, self%nsppol, self%ndat) 6054 6055 !Local variables------------------------------- 6056 integer :: spin, idat, band 6057 !integer :: ir1, ir2, ir3, pr1, pr2, pr3 6058 real(dp) :: val !, vv(8) 6059 real(dp) :: kwrap(3), shift(3) 6060 6061 ! ********************************************************************* 6062 6063 call wrap2_zero_one(kpt, kwrap, shift) 6064 !write(std_out, *)"kwrap:", kwrap 6065 6066 ! ir1,ir2,ir3 = bottom left neighbor 6067 ! pr1,pr2,pr3 = top right neighbor 6068 !call interpol3d_indices(kwrap, self%nkx, self%nky, self%nkz, ir1, ir2, ir3, pr1, pr2, pr3) 6069 6070 do idat=1,self%ndat 6071 do spin=1,self%nsppol 6072 do band=1,self%bsize 6073 val = interpol3d_0d(kwrap, self%nkx, self%nky, self%nkz, self%data_uk_bsd(:,:,:,band, spin, idat)) 6074 6075 !if (val <= zero) then 6076 ! vv(1) = self%data_uk_bsd(ir1, ir2, ir3, band, spin, idat) 6077 ! vv(2) = self%data_uk_bsd(pr1, ir2, ir3, band, spin, idat) 6078 ! vv(3) = self%data_uk_bsd(ir1, pr2, ir3, band, spin, idat) 6079 ! vv(4) = self%data_uk_bsd(ir1, ir2, pr3, band, spin, idat) 6080 ! vv(5) = self%data_uk_bsd(pr1, pr2, ir3, band, spin, idat) 6081 ! vv(6) = self%data_uk_bsd(ir1, pr2, pr3, band, spin, idat) 6082 ! vv(7) = self%data_uk_bsd(pr1, ir2, pr3, band, spin, idat) 6083 ! vv(8) = self%data_uk_bsd(pr1, pr2, pr3, band, spin, idat) 6084 ! val = maxval(vv) 6085 !end if 6086 6087 vals_bsd(band, spin, idat) = val 6088 end do 6089 end do 6090 end do 6091 6092 end subroutine klinterp_eval_bsd
m_ebands/klinterp_free [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
klinterp_free
FUNCTION
Free dynamic memory.
INPUTS
OUTPUT
SOURCE
6020 subroutine klinterp_free(self) 6021 6022 !Arguments ------------------------------------ 6023 !scalars 6024 class(klinterp_t),intent(inout) :: self 6025 6026 ! ********************************************************************* 6027 6028 ABI_SFREE(self%data_uk_bsd) 6029 6030 end subroutine klinterp_free
m_ebands/klinterp_new [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
klinterp_new
FUNCTION
INPUTS
OUTPUT
SOURCE
5901 type(klinterp_t) function klinterp_new(cryst, kptrlatt, nshiftk, shiftk, kptopt, kibz, & 5902 bsize, nkibz, nsppol, ndat, values_bksd, comm) result(new) 5903 5904 !Arguments ------------------------------------ 5905 !scalars 5906 type(crystal_t),intent(in) :: cryst 5907 integer,intent(in) :: nshiftk, kptopt, bsize, nkibz, nsppol, ndat, comm 5908 !arrays 5909 integer,intent(in) :: kptrlatt(3,3) 5910 real(dp),intent(in) :: kibz(3, nkibz), shiftk(3,nshiftk), values_bksd(bsize, nkibz, nsppol, ndat) 5911 5912 !Local variables------------------------------- 5913 !scalars 5914 integer,parameter :: sppoldbl1 = 1 5915 integer :: ierr, nkfull, ikf, ik_ibz, timrev, ix, iy, iz, nkx, nky, nkz !spin, band, idat 5916 real(dp) :: dksqmax 5917 character(len=500) :: msg 5918 !arrays 5919 integer,allocatable :: bz2ibz(:,:) 5920 real(dp) :: kpt(3) 5921 real(dp),allocatable :: kfull(:,:) 5922 5923 ! ********************************************************************* 5924 5925 ! Check input parameters 5926 ierr = 0 5927 if (nkibz == 1) then 5928 ABI_ERROR_NOSTOP("Cannot interpolate with a single k-point", ierr) 5929 end if 5930 if (.not. isdiagmat(kptrlatt)) then 5931 ABI_ERROR_NOSTOP('kptrlatt is not diagonal. Multiple shifts are not allowed', ierr) 5932 end if 5933 if (nshiftk /= 1) then 5934 ABI_ERROR_NOSTOP('Multiple shifts not allowed', ierr) 5935 end if 5936 if (any(abs(shiftk(:, 1)) > tol8)) then 5937 ABI_ERROR_NOSTOP("shifted k-mesh not implented", ierr) 5938 end if 5939 5940 if (ierr /= 0) then 5941 ABI_ERROR("Linear interpolation cannot be performed. See messages above.") 5942 end if 5943 5944 nkx = kptrlatt(1, 1) 5945 nky = kptrlatt(2, 2) 5946 nkz = kptrlatt(3, 3) 5947 5948 new%nkx = nkx; new%nky = nky; new%nkz = nkz 5949 new%bsize = bsize; new%nsppol = nsppol; new%ndat = ndat 5950 5951 ! Build list of k-points in the conventional unit cell. 5952 ! (x,y,z) ordered as required by interpolation routine 5953 nkfull = nkx * nky * nkz 5954 ABI_MALLOC(kfull, (3, nkfull)) 5955 ikf = 0 5956 do iz=1,nkz 5957 kpt(3) = (iz - 1 + shiftk(3, 1)) / nkz 5958 do iy=1,nky 5959 kpt(2) = (iy - 1 + shiftk(2, 1)) / nky 5960 do ix=1,nkx 5961 kpt(1) = (ix - 1 + shiftk(1, 1)) / nkx 5962 ikf = ikf + 1 5963 kfull(:, ikf) = kpt 5964 end do 5965 end do 5966 end do 5967 5968 ! Build mapping kfull --> IBZ 5969 timrev = kpts_timrev_from_kptopt(kptopt) 5970 ABI_MALLOC(bz2ibz, (nkfull*sppoldbl1, 6)) 5971 5972 call listkk(dksqmax, cryst%gmet, bz2ibz, kibz, kfull, nkibz, nkfull, cryst%nsym,& 5973 sppoldbl1, cryst%symafm, cryst%symrec, timrev, comm, use_symrec=.True.) 5974 5975 ABI_FREE(kfull) 5976 5977 if (dksqmax > tol12) then 5978 write(msg, '(3a,es16.6,4a)' )& 5979 'At least one of the k points could not be generated from a symmetrical one.',ch10,& 5980 'dksqmax: ',dksqmax,ch10,& 5981 'Action: check k-point input variables',ch10,& 5982 ' e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.' 5983 ABI_ERROR(msg) 5984 end if 5985 5986 ABI_CALLOC(new%data_uk_bsd, (nkx, nky, nkz, bsize, nsppol, ndat)) 5987 5988 ! Build array in the full BZ to prepare call to interpol3d_0d. 5989 ikf = 0 5990 do iz=1,nkz 5991 do iy=1,nky 5992 do ix=1,nkx 5993 ikf = ikf + 1 5994 ik_ibz = bz2ibz(ikf, 1) 5995 new%data_uk_bsd(ix, iy, iz, 1:bsize, 1:nsppol, 1:ndat) = values_bksd(1:bsize, ik_ibz, 1:nsppol, 1:ndat) 5996 end do 5997 end do 5998 end do 5999 6000 ABI_FREE(bz2ibz) 6001 6002 end function klinterp_new
m_ebands/klinterp_t [ Types ]
[ Top ] [ m_ebands ] [ Types ]
NAME
klinterp_t
FUNCTION
Linear interpolation of eigenvalue-like quantities (scalars with the same symmetry as the KS eigenvalues) Used, for instance, to interpolate electron or phonon lifetimes.
SOURCE
341 type,public :: klinterp_t 342 343 integer :: bsize, nsppol, ndat 344 ! Max number of bands, number of independent spin polarization, size of "extra" dimension. 345 346 integer :: nkx, nky, nkz 347 ! Number of divisions of the grid enclosing the first unit cell 348 349 real(dp),allocatable :: data_uk_bsd(:,:,:,:,:,:) 350 ! (nkx*nky*nkz, mband, nsppol, ndat) 351 352 contains 353 354 procedure :: free => klinterp_free 355 ! Free dynamic memory 356 357 procedure :: eval_bsd => klinterp_eval_bsd 358 ! Interpolate values at an arbitrary k-point. 359 360 end type klinterp_t 361 362 public :: klinterp_new ! Build interpolator. 363 364 365 !---------------------------------------------------------------------- 366 367 CONTAINS !=====================================================================================
m_ebands/pack_eneocc [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
pack_eneocc
FUNCTION
Helper function to do a reshape of (energies|occupancies|derivate of occupancies) initially stored in a 3D arrays returning a vector.
INPUTS
nkpt=number of k-points nsppol=number of spin polarizations mband=Max number of bands over k-points (just to dimension the output) nbands(nkpt*nsppol)=Number of bands at eack k and spin bantot=Total number of bands array3d(mband,nkpt,nsppol)=Arrays containing the values to reshape.
OUTPUT
vect(bantot)=The input values stored in vector mode. Only the values really considered at each k-point and spin are copied.
SOURCE
1323 subroutine pack_eneocc(nkpt, nsppol, mband, nband, bantot, array3d, vect) 1324 1325 !Arguments ------------------------------------ 1326 !scalars 1327 integer,intent(in) :: nkpt,nsppol,mband,bantot 1328 !arrays 1329 integer,intent(in) :: nband(nkpt*nsppol) 1330 real(dp),intent(in) :: array3d(mband,nkpt,nsppol) 1331 real(dp),intent(out) :: vect(bantot) 1332 1333 !Local variables------------------------------- 1334 integer :: spin,ikpt,band,idx 1335 1336 ! ************************************************************************* 1337 1338 vect(:)=zero 1339 idx=0 1340 do spin=1,nsppol 1341 do ikpt=1,nkpt 1342 do band=1,nband(ikpt+(spin-1)*nkpt) 1343 idx=idx+1 1344 vect(idx)=array3d(band,ikpt,spin) 1345 end do 1346 end do 1347 end do 1348 1349 end subroutine pack_eneocc
m_ebands/put_eneocc_vect [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
put_eneocc_vect
FUNCTION
Update the energies or the occupations stored in a ebands_t structure. The input values are stored in a vector according to the abinit convention In the data type, on the contrary, we use 3D arrays (mband,nkpt,nsspol) which are much easier to use inside loops.
INPUTS
vect(ebands%bantot)=The new values to be stored in the structure. arr_name=The name of the quantity to be saved (CASE insensitive). Allowed values are == "eig" == For the eigenvalues. == "occ" == For the occupation numbers. == "doccde" == For the derivative of the occupancies wrt the energy.
OUTPUT
See SIDE EFFECTS
SIDE EFFECTS
ebands<ebands_t>=The object with updated values depending on the value of arr_name
SOURCE
1430 subroutine put_eneocc_vect(ebands,arr_name,vect) 1431 1432 !Arguments ------------------------------------ 1433 !scalars 1434 character(len=*),intent(in) :: arr_name 1435 class(ebands_t),intent(inout) :: ebands 1436 real(dp),intent(in) :: vect(:) 1437 1438 !Local variables------------------------------- 1439 integer :: nkpt,nsppol,mband,bantot 1440 real(dp) :: val 1441 ! ************************************************************************* 1442 1443 mband =ebands%mband; bantot=ebands%bantot; nkpt =ebands%nkpt; nsppol=ebands%nsppol 1444 1445 select case (tolower(arr_name)) 1446 case ('occ') 1447 call unpack_eneocc(nkpt,nsppol,mband,ebands%nband,vect,ebands%occ, val=zero) 1448 case ('eig') 1449 ! DFPT routines call ebands_init with the wrong bantot. Using maxval(vect) causes SEGFAULT 1450 ! so I have to recompute the correct bantot here 1451 !ABI_CHECK(sum(ebands%nband) == ebands%bantot, "bantot and nband are incosistent") 1452 val = maxval(vect(1:sum(ebands%nband))) 1453 call unpack_eneocc(nkpt,nsppol,mband,ebands%nband,vect,ebands%eig, val=val) 1454 case ('doccde') 1455 call unpack_eneocc(nkpt,nsppol,mband,ebands%nband,vect,ebands%doccde, val=zero) 1456 case default 1457 ABI_BUG(sjoin('Wrong arr_name= ', arr_name)) 1458 end select 1459 1460 end subroutine put_eneocc_vect
m_ebands/unpack_eneocc [ Functions ]
[ Top ] [ m_ebands ] [ Functions ]
NAME
unpack_eneocc
FUNCTION
Helper function to do a reshape of (energies|occupancies|derivate of occupancies) initially stored in a vector. Return a 3D array index by (band,ikpt,spin)
INPUTS
nkpt=number of k-points nsppol=number of spin polarizations mband=Max number of bands over k-points (just to dimension the output) nbands(nkpt*nsppol)=Number of bands at eack k and spin vect(:)=The input values to reshape [val]=Optional value used to initialize the array.
OUTPUT
array3d(mband,nkpt,nsppol)=Arrays containing the values of vect. Note that the first dimension is usually larger than the number of bands really used for a particular k-point and spin.
SOURCE
1264 subroutine unpack_eneocc(nkpt,nsppol,mband,nband,vect,array3d,val) 1265 1266 !Arguments ------------------------------------ 1267 !scalars 1268 integer,intent(in) :: nkpt,nsppol,mband 1269 real(dp),optional,intent(in) :: val 1270 !arrays 1271 integer,intent(in) :: nband(nkpt*nsppol) 1272 real(dp),intent(in) :: vect(:) 1273 real(dp),intent(out) :: array3d(mband,nkpt,nsppol) 1274 1275 !Local variables------------------------------- 1276 integer :: spin,ikpt,band,idx 1277 ! ************************************************************************* 1278 1279 if (present(val)) then 1280 array3d = val 1281 else 1282 array3d = huge(one) 1283 end if 1284 1285 idx=0 1286 ! elements in vect are packed in the first positions. 1287 do spin=1,nsppol 1288 do ikpt=1,nkpt 1289 do band=1,nband(ikpt + (spin-1)*nkpt) 1290 idx = idx + 1 1291 array3d(band, ikpt, spin) = vect(idx) 1292 end do 1293 end do 1294 end do 1295 1296 end subroutine unpack_eneocc