1!==============================================================================
2!
3! Module:
4!
5! (1) read_wannier_m        Originally By DAS  April 2012
6!
7!     Reads binary checkpoint file "seedname.chk" produced by Wannier90 v1.2
8!     to obtain the coefficients relating the Wannier to Bloch functions.
9!     Reference: wannier90-1.2/src/parameters.F90 routine param_write_chkpt
10!
11!==============================================================================
12
13#include "f_defs.h"
14
15module read_wannier_m
16
17  use global_m
18  implicit none
19
20  private
21
22  public :: read_wannier
23
24contains
25
26  subroutine read_wannier(seedname, crys, kg, iwann, nvband, coefficients)
27    character(len=*), intent(in) :: seedname !< read file seedname.chk
28    type(crystal), intent(in) :: crys !< for checking lattice vectors
29    type(grid), intent(in) :: kg !< for assigning k-points
30    integer, intent(in) :: iwann !< index of Wannier function to read
31    integer, intent(in) :: nvband !< number of valence bands in plotxct calculation
32    complex(DPC), intent(out) :: coefficients(:,:) !< u_matrix relating Bloch to Wannier. (band, k-point)
33
34    integer, parameter :: chk_unit = 73
35    integer :: ib,iw,ii,jj,ik,ik2,num_bands,num_kpts,num_wann,ikmatch
36    character :: filename*50, header*33, chkpt1*20
37    logical :: have_disentangled
38    real(DP) :: real_lattice(3,3), recip_lattice(3,3)
39    real(DP), allocatable :: kpt_latt(:,:)
40    complex(DPC), allocatable :: u_matrix(:,:,:)
41
42    PUSH_SUB(read_wannier)
43
44    filename = trim(seedname)//'.chk'
45    write(6,'(a)') "Reading file " // trim(filename)
46    call open_file(unit=chk_unit,file=filename,form='unformatted',status='old')
47
48    if(iwann <= 0) then
49      call die("iwann cannot be negative", only_root_writes = .true.)
50    endif
51
52    read(chk_unit) header                                       ! Date and time
53    write(6,'(a)') header
54
55    read(chk_unit) num_bands                                    ! Number of bands
56    if(num_bands < nvband) then
57      if(peinf%inode == 0) then
58        write(0,'(a,9f11.6)') 'read ', num_bands, ' bands but need at least ', nvband, ' = number of valence bands'
59      endif
60      call die("not enough bands in " // trim(filename), only_root_writes = .true.)
61    endif
62
63    read(chk_unit) ! num_exclude_bands (don`t care)             ! Number of excluded bands
64    read(chk_unit) ! (exclude_bands(ib),ib=1,num_exclude_bands) ! Excluded bands
65
66    read(chk_unit) ((real_lattice(ii,jj),ii=1,3),jj=1,3)        ! Real lattice (a.u.)
67    if(any(abs(real_lattice(:,:) - crys%avec(:,:) * crys%alat) > TOL_Small)) then
68      if(peinf%inode == 0) then
69        write(0,'(a,9f11.6)') 'read lattice vectors ', real_lattice(:,:)
70        write(0,'(a,9f11.6)') 'should be ', crys%avec(:,:) * crys%alat
71      endif
72      call die("incorrect real-space lattice vectors", only_root_writes = .true.)
73    endif
74
75    read(chk_unit) ((recip_lattice(ii,jj),ii=1,3),jj=1,3)       ! Reciprocal lattice (a.u.)
76    if(any(abs(real_lattice(:,:) - crys%bvec(:,:) * crys%alat) > TOL_Small)) then
77      if(peinf%inode == 0) then
78        write(0,'(a,9f11.6)') 'read lattice vectors ', recip_lattice(:,:)
79        write(0,'(a,9f11.6)') 'should be ', crys%bvec(:,:) * crys%blat
80      endif
81      call die("incorrect reciprocal lattice vectors", only_root_writes = .true.)
82    endif
83
84    read(chk_unit) num_kpts                                     ! Number of k-points
85    if(num_kpts < kg%nf) then
86      if(peinf%inode == 0) write(0,'(a,i6,a,i6)') 'Need ', kg%nf, 'kpoints but file has only ', num_kpts
87      call die("Not enough kpoints in " // trim(filename), only_root_writes = .true.)
88    endif
89
90    read(chk_unit) !(mp_grid(ii),ii=1,3) (don`t care)           ! M-P grid
91    SAFE_ALLOCATE(kpt_latt, (3, num_kpts))
92    read(chk_unit) ((kpt_latt(ii,ik),ii=1,3),ik=1,num_kpts)     ! K-points in lattice vectors
93    read(chk_unit) !nntot (don`t care)                          ! Number of nearest k-point neighbours
94
95    read(chk_unit) num_wann                                     ! Number of wannier functions
96    if(iwann > num_wann) then
97      if(peinf%inode == 0) write(0,'(2a)') 'Requested Wannier function ', iwann, ' > number of Wannier functions ', num_wann
98      call die("iwann out of bounds", only_root_writes = .true.)
99    endif
100    read(chk_unit) chkpt1                                       ! Position of checkpoint
101    if(trim(chkpt1) /= 'post_wann') then
102      if(peinf%inode == 0) write(0,'(2a)') 'chkpt1 = ', trim(chkpt1)
103      call die("Wannier90 checkpoint position must be post_wann.", only_root_writes = .true.)
104    endif
105    read(chk_unit) have_disentangled                            ! Whether a disentanglement has been performed
106    if (have_disentangled) then
107      call die("disentanglement not implemented")
108!      read(chk_unit) omega_invariant     ! Omega invariant
109! lwindow, ndimwin and U_matrix_opt
110!      read(chk_unit) ((lwindow(i,nkp),i=1,num_bands),nkp=1,num_kpts)
111!      read(chk_unit) (ndimwin(nkp),nkp=1,num_kpts)
112!      read(chk_unit) (((u_matrix_opt(i,j,nkp),i=1,num_bands),j=1,num_wann),nkp=1,num_kpts)
113    endif
114    SAFE_ALLOCATE(u_matrix, (num_wann, num_wann, num_kpts))
115    read(chk_unit) (((u_matrix(ib,iw,ik),ib=1,num_wann),iw=1,num_wann),ik=1,num_kpts) ! U_matrix
116!    read(chk_unit) ((((m_matrix(i,j,k,l),i=1,num_wann),j=1,num_wann),k=1,nntot),l=1,num_kpts) ! M_matrix
117!    read(chk_unit) ((wannier_centres(i,j),i=1,3),j=1,num_wann)
118!    read(chk_unit) (wannier_spreads(ib),i=1,num_wann)
119    call close_file(chk_unit)
120
121    ! relate wannier90 kpts to those in plotxct
122    do ik = 1, kg%nf
123      ikmatch = 0
124      do ik2 = 1, num_kpts
125        if(all(abs(kg%f(1:3, ik) - kpt_latt(1:3, ik2)) < TOL_Small)) then
126          ikmatch = ik2
127          exit
128        endif
129      enddo
130      if(ikmatch > 0) then
131        do ib = 1, nvband
132          coefficients(ib, ik) = u_matrix(nvband - ib + 1, iwann, ikmatch)
133        enddo
134      else
135        if(peinf%inode == 0) write(0,*) 'no match for kpoint ', kg%f(1:3, ik)
136        call die("need kpoint not present in " // trim(filename), only_root_writes = .true.)
137      endif
138    enddo
139
140    SAFE_DEALLOCATE(u_matrix)
141    SAFE_DEALLOCATE(kpt_latt)
142
143    POP_SUB(read_wannier)
144    return
145  end subroutine read_wannier
146
147end module read_wannier_m
148