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