1!----------------------------------------------------------------------- 2! 3! module fullbz_m 4! 5! Routines: fullbz 6! 7! If wigner_seitz = .true. (for BSE, PlotXct, NonLinearOptics) 8! Uses a Wigner-Seitz construction to define the Brillouin zone. 9! If wigner_seitz = .false. (for Epsilon, Sigma) 10! Uses the usual "box" BZ. 11! 12! input: crys,syms type 13! ntran (number of symmetry operations) 14! gr%nr 15! gr%rk 16! 17! output: gr type (except gr%nr, gr%rk) 18! 19!----------------------------------------------------------------------- 20 21#include "f_defs.h" 22 23module fullbz_m 24 25 use global_m 26 use misc_m 27 implicit none 28 29 private 30 31 public :: fullbz, dealloc_grid 32 33contains 34 35 !> Set up the mapping between the irreducible and the full full Brillouin zone 36 !! i.e. the k-point correspondence and the symmetry operations. 37 !! Everything is stored into the grid object. 38 subroutine fullbz(crys,syms,gr,ntran,skip_checkbz,wigner_seitz,paranoid,do_nothing,nfix) 39 type (crystal), intent(in) :: crys 40 type (symmetry), intent(in) :: syms 41 type (grid), intent(inout) :: gr 42 integer, intent(in) :: ntran !< number of sym ops, typically syms%ntrans 43 logical, intent(out) :: skip_checkbz 44 logical, intent(in) :: wigner_seitz !< do a Wigner-Seitz construction 45 logical, intent(in) :: paranoid !< perform paranoia check: no k-pts differ by G-vector 46 logical, optional, intent(in) :: do_nothing !< initialize without performing any checks. 47 integer, optional, intent(in) :: nfix !< don`t unfold the first nfix points 48 49 integer :: ii,jj,kk,it,ir,if,i1,i2,i3,gpt(3),iostat_c,ntran_ 50 real(DP) :: tmpf(3),tmpfm(3),length,lmin,fq(3) 51 real(DP), allocatable :: fk(:,:),kg0(:,:) 52 integer, allocatable :: indr(:),itran(:) 53 logical :: found 54! 55! i r,nr k-points in irr bz 56! o fk,gr%nf k-points in full bz 57! o sz radius of a spherical subzone equivalent to 58! one point in the set fk 59! 60! Loop over ir-points 61! 62 PUSH_SUB(fullbz) 63 64 if(present(do_nothing)) then 65 if(do_nothing) then 66 if(wigner_seitz .or. paranoid) call die("bug: cannot call fullbz with do_nothing and wigner_seitz or paranoid.") 67 gr%nf=gr%nr 68 gr%sz=2.0d0*PI_D*(3.0d0/(4.0d0*PI_D*gr%nf*crys%celvol))**(1.0d0/3.0d0) 69 SAFE_ALLOCATE(gr%kg0,(3,gr%nf)) 70 SAFE_ALLOCATE(gr%f,(3,gr%nf)) 71 SAFE_ALLOCATE(gr%itran,(gr%nf)) 72 SAFE_ALLOCATE(gr%indr,(gr%nf)) 73 gr%kg0(:,:)=0 74 gr%itran(:)=1 75 gr%f(:,:)=gr%r(:,:) 76 do ii=1,gr%nf 77 gr%indr(ii)=ii 78 enddo 79 skip_checkbz = .true. 80 POP_SUB(fullbz) 81 return 82 endif 83 endif 84 85 SAFE_ALLOCATE(fk, (3,gr%nr*ntran)) 86 SAFE_ALLOCATE(indr, (gr%nr*ntran)) 87 SAFE_ALLOCATE(itran, (gr%nr*ntran)) 88 89 call open_file(unit=21,file='fullbz.dat',status='old',iostat=iostat_c) 90 skip_checkbz = (iostat_c==0) 91 if (skip_checkbz) then 92 write(6,'(3x,a)') 'Reading the full Brillouin Zone from fullbz.dat' 93 write(6,'(3x,a)') 'Will not be checking if the points there form a full zone' 94 read(21,*) gr%nf 95 do ii=1,gr%nf 96 read(21,*) fk(1:3,ii),itran(ii),indr(ii) 97 enddo 98 call close_file(21) 99 ! This is just in case we want to do wigner_seitz or paranoid 100 else 101 gr%nf=0 102 103 do ir=1,gr%nr 104! 105! Loop over transformations 106! 107 ntran_ = ntran 108 if (present(nfix)) then 109 if (ir<=nfix) ntran_ = 1 110 endif 111 do it=1,ntran_ 112! 113! Rotate gr%r and put into tmpf 114! 115 tmpf(:) = matmul(dble(syms%mtrx(:,:,it)),gr%r(:,ir)) 116 117 call k_range(tmpf, gpt, TOL_Small) 118! 119! Compare to other points in full zone 120! 121 found = .false. 122 do if=1,gr%nf 123 if(all(abs(tmpf(1:3)-fk(1:3,if)).lt.TOL_Small)) then 124 found = .true. 125 exit 126 endif 127 enddo 128 if(found) cycle ! skip adding it 129! 130! Store new kpoint in fbz 131! 132 gr%nf=gr%nf+1 133 if (gr%nf > gr%nr * ntran) call die('fullbz internal error') 134 fk(1:3,gr%nf)=tmpf(1:3) 135! 136! Store index of rotation itran and corresponding IBZ point 137! 138 itran(gr%nf)=it 139 indr(gr%nf)=ir 140 141 enddo !end loop over symmetries 142 enddo !end loop over the q-points from epsmat and eps0mat 143 endif ! whether reading from fullbz.dat 144 145! 146! SIB: Now that we have the full BZ with components [0,1), we 147! will move each point into the Wigner-Seitz cell by adding G-vectors 148! to it until we get the shortest length vector we can. Then we find 149! the appropriate Umklapp vector as well (kg0). 150! 151 152 if(wigner_seitz) then 153 SAFE_ALLOCATE(kg0, (3,gr%nr*ntran)) 154 155 do ii=1,gr%nf 156 tmpf(:) = fk(:,ii) 157 lmin = 1.0d10 158 do i1=-ncell+1,ncell 159 do i2=-ncell+1,ncell 160 do i3=-ncell+1,ncell 161 fq(1) = tmpf(1) - i1 162 fq(2) = tmpf(2) - i2 163 fq(3) = tmpf(3) - i3 164 length = DOT_PRODUCT(fq,MATMUL(crys%bdot,fq)) 165 if (length.lt.lmin) then 166 lmin = length 167 tmpfm(:) = fq(:) 168 endif 169 enddo 170 enddo 171 enddo 172 fk(:,ii) = tmpfm 173! 174! SIB: Find Umklapp (kg0) a.k.a. translation 175! 176 tmpf = MATMUL(dble(syms%mtrx(:,:,itran(ii))),gr%r(:,indr(ii))) 177 tmpf(:) = fk(:,ii) - tmpf 178 do jj=1,3 179 if (tmpf(jj).ge.0.0) kg0(jj,ii)=tmpf(jj)+TOL_Small 180 if (tmpf(jj).lt.0.0) kg0(jj,ii)=tmpf(jj)-TOL_Small 181 enddo 182 183 enddo !over full zone (ii) 184 185 SAFE_ALLOCATE(gr%kg0, (3,gr%nf)) 186 gr%kg0(1:3,1:gr%nf)=kg0(1:3,1:gr%nf) 187 SAFE_DEALLOCATE(kg0) 188 189 else 190 nullify(gr%kg0) 191 endif 192 193!------------------------------------------------------------------------ 194! SIB: Paranoia check. We will compare all pairs of points in the full 195! zone to each other in order to see if any of them differ by a G-vector. 196! If they do, we are in big trouble and stop! 197 198 199 if(paranoid) then 200 ii_loop: do ii=1,gr%nf 201 do jj=1,ii-1 202 tmpf = abs(fk(:,ii)-fk(:,jj)) 203! after this loop, tmpf contains how far each component of fii-fjj is 204! from the closest integer to it. 205 do kk=1,3 206 it = tmpf(kk) 207 tmpf(kk) = tmpf(kk)-it 208 if (tmpf(kk).ge.0.5) tmpf(kk)=1.0-tmpf(kk) 209 enddo 210! if fkii-fkjj is very close to an integer vector (G-vector), trouble! 211 if (sum(abs(tmpf)).le.TOL_Small) then 212 write(0,123) ii,jj,fk(:,ii)-fk(:,jj) 213123 format('equiv kpts',i4,' and ',i4,' diff=',3f10.5) 214 call die('equivalent points found in the full BZ') 215 endif 216 enddo 217 enddo ii_loop 218 endif 219 220!------------------------------------------------------------------------ 221! SIB: Now store all gathered information where it belongs 222 223 SAFE_ALLOCATE(gr%f, (3,gr%nf)) 224 SAFE_ALLOCATE(gr%indr, (gr%nf)) 225 SAFE_ALLOCATE(gr%itran, (gr%nf)) 226 gr%f(1:3,1:gr%nf)=fk(1:3,1:gr%nf) 227 gr%indr(1:gr%nf)=indr(1:gr%nf) 228 gr%itran(1:gr%nf)=itran(1:gr%nf) 229 SAFE_DEALLOCATE(fk) 230 SAFE_DEALLOCATE(indr) 231 SAFE_DEALLOCATE(itran) 232 233! Compute radius of spherical subzone 234! assuming bz filled with gr%nf spheres 235 236 gr%sz=2.0d0*PI_D*(3.0d0/(4.0d0*PI_D*gr%nf*crys%celvol))**(1.0d0/3.0d0) 237 238 if (ntran.eq.1 .and. gr%nf.ne.gr%nr .and. paranoid) then 239 write(0,*) gr%nf, ' and ', gr%nr 240 call die('fullbz paranoia check failed') 241 endif 242 243 POP_SUB(fullbz) 244 245 return 246 end subroutine fullbz 247 248 !----------------------------------------------------------------------- 249 subroutine dealloc_grid(gr) 250 type(grid), intent(inout) :: gr 251 252 PUSH_SUB(dealloc_grid) 253 254 SAFE_DEALLOCATE_P(gr%r) 255 SAFE_DEALLOCATE_P(gr%f) 256 SAFE_DEALLOCATE_P(gr%indr) 257 SAFE_DEALLOCATE_P(gr%itran) 258 SAFE_DEALLOCATE_P(gr%kg0) 259 260 POP_SUB(dealloc_grid) 261 return 262 end subroutine dealloc_grid 263 264end module fullbz_m 265