1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8module m_siesta2wannier90 9 10! OUTPUT: 11! File called seedname.mmn, where the overlap matrices are written in 12! the format required by Wannier90 13! 14! File called seedname.eig, where the eigenvalues of the Hamiltonian 15! at the required k-points are written according to the format required 16! by Wannier90. 17 18 use precision, only : dp ! Real double precision type 19 use parallel, only : IOnode ! Input/output node 20 use parallel, only : Node ! This process node 21 use parallel, only : Nodes ! Total number of processor nodes 22 use sys, only : die ! Termination routine 23 use files, only : label_length ! Number of characters in slabel 24 use siesta_options, only: w90_write_mmn ! Write the Mmn matrix for the 25 ! interface with Wannier 26 use siesta_options, only: w90_write_amn ! Write the Amn matrix for the 27 ! interface with Wannier 28 use siesta_options, only: w90_write_eig ! Write the eigenvalues for the 29 ! interface with Wannier 30 use siesta_options, only: w90_write_unk ! Write the unks for the 31 ! interface with Wannier 32 use TrialOrbitalClass 33 34! 35! Variables related with the atomic structure of the system 36! 37 38 real(dp) :: latvec(3,3) ! Real lattice vectors. 39 ! Cartesian coordinates. 40 ! Readed in Angstroms and transformed to Bohr 41 ! internally 42 ! First index: component 43 ! Second index: vector 44 ! This is consistent with the unit cell read 45 ! in Siesta, but it has changed with respect 46 ! the first version implemented by R. Korytar 47 real(dp) :: reclatvec(3,3) ! Reciprocal lattice vectors. 48 ! Cartesian coordinates. 49 ! Readed in Angstroms^-1 and transformed 50 ! to Bohr^-1 internally 51 ! First index: component 52 ! Second index: vector 53 ! This is consistent with the reciprocal 54 ! lattice read in Siesta 55 ! in Siesta, but it has changed with respect 56 ! the first version implemented by R. Korytar 57! 58! Variables related with the k-point list for which the overlap 59! matrices Mmn between a k-point and its neighbor will be computed 60! 61 integer :: numkpoints ! Total number of k-points 62 ! for which the overlap of the 63 ! periodic part of the wavefunct 64 ! with a neighbour k-point will 65 ! be computed 66 real(dp), pointer :: kpointsfrac(:,:) 67 ! List of k points relative 68 ! to the reciprocal lattice vectors. 69 ! First index: component 70 ! Second index: k-point index in the list 71 real(dp), pointer :: bvectorsfrac(:,:) 72 ! The vectors b that connect each mesh-point k 73 ! to its nearest neighbours 74! 75! Variables related with the neighbours of the k-points 76! 77 integer :: nncount 78 ! The number of nearest neighbours belonging to 79 ! each k-point of the Monkhorst-Pack mesh 80 integer, pointer :: nnlist(:,:) 81 ! nnlist(ikp,inn) is the index of the 82 ! inn-neighbour of ikp-point 83 ! in the Monkhorst-Pack grid folded to the 84 ! first Brillouin zone 85 integer, pointer :: nnfolding(:,:,:) 86 ! nnfolding(i,ikp,inn) is the i-component 87 ! of the reciprocal lattice vector 88 ! (in reduced units) that brings 89 ! the inn-neighbour specified in nnlist 90 ! (which is in the first BZ) to the 91 ! actual \vec{k} + \vec{b} that we need. 92 ! In reciprocal lattice units. 93! 94! Variables related with the projections with trial functions, 95! initial approximations to the MLWF 96! 97 integer :: numproj ! Total number of projection centers, 98 ! equal to the number of MLWF 99 type(trialorbital), target, allocatable :: projections(:) 100! 101! Variables related with the number of bands considered for wannierization 102! 103 integer :: numbands(2) ! Number of bands for wannierization 104 ! before excluding bands 105 integer :: numexcluded 106 ! Number of bands to exclude from the calculation 107 ! of the overlap and projection matrices. 108 ! This variable is read from the .nnkp file 109 110 integer :: numincbands(2) 111 ! Number of included bands in the calc. 112 ! of the overlap and projection matrices. 113 ! after excluding the bands 114 integer :: nincbands_loc 115 ! Number of included bands in the calc. 116 ! of the overlap and projection matrices. 117 ! in the local Node 118 integer :: blocksizeincbands ! Maximum number of bands 119 ! considered for wannierization per node 120 121 integer, pointer :: excludedbands(:) 122 ! Bands to be excluded 123 ! This variable is read from the .nnkp file 124 logical, pointer :: isexcluded(:) ! Masks excluded bands 125 integer, pointer :: isincluded(:) ! Masks included bands 126 127! 128! Variables related with the coefficients of the wavefunctions and 129! eigenvalues at the Wannier90 k-point mesh 130! 131 complex(dp), pointer :: coeffs(:,:,:) => null() ! Coefficients of the wavefunctions. 132 ! First index: orbital 133 ! Second index: band 134 ! Third index: k-point 135 real(dp), pointer :: eo(:,:) => null() ! Eigenvalues of the Hamiltonian 136 ! at the numkpoints introduced in 137 ! kpointsfrac 138 ! First index: band index 139 ! Second index: k-point index 140 141 142! 143! Output matrices 144! 145 146 complex(dp), pointer :: Mmnkb(:,:,:,:) => null() ! Matrix of the overlaps of 147 ! periodic parts of Bloch waves. 148 ! <u_{ik}|u_{jk+b}> 149 ! The first two indices refer to 150 ! the number of occupied bands 151 ! (indices m and n in standard 152 ! notation, see for instance, 153 ! Eq. (27) of the paper by 154 ! Marzari et al., RMP 84, 1419 (2012) 155 ! The third index refer to the kpoint 156 ! The fourth index refer to the neig 157 complex(dp), pointer :: Amnmat(:,:,:) => null() ! Projections of a trial function 158 ! with a Bloch orbital 159 ! <\psi_{m k}|g_n> 160 161! 162! Variables related with the input/output files 163! 164 character(label_length+3) :: seedname ! Name of the file where the Wannier90 165 ! code, when used as a postprocessing 166 ! tool, reads or dumps the 167 ! information. 168 169CONTAINS 170 171subroutine siesta2wannier90 172 173 use m_spin, only: nspin ! Number of spin components 174 use files, only: slabel ! Short system label, 175 ! used to generate file names 176 use m_digest_nnkp, only: read_nnkp ! Subroutine that reads the .nnkp file 177 use m_digest_nnkp, only: chosing_b_vectors ! Subroutine that computes the b 178 ! vectors that connect each mesh 179 ! k-point to its nearest neighbours. 180 use m_digest_nnkp, only: set_excluded_bands ! Subroutine that chooses the 181 ! bands that are excluded from 182 ! the wannierization procedure 183 use m_digest_nnkp, only: number_bands_wannier ! Subroutine that computes the 184 ! number of bands for wannierization 185 186 implicit none 187 188 integer :: ispin ! Spin counter 189 integer :: numbandswan(2) 190 ! Number of bands for wannierization 191 192 external :: timer, mmn 193 194 call timer("siesta2wannier90",1) 195 196 do ispin = 1, nspin 197! Append _up or _dn to the seedname if spin polarized 198 seedname = trim(getFileNameRoot(ispin,nspin,slabel)) 199! A priori, the .nnkp file generated by Wannier90 used as a 200! postprocessing tool is independent of spin. 201! However, they provide examples (example08, bulk Fe) where 202! two different .nnkp files are generated, 203! one for each component of the spin. 204! This inspired us to read the .nnkp file for each component of the spin 205! (note that the seedname will be different in every case). 206! The information is read only by the master node, 207! and broadcast to the rest of the nodes later. 208 if (IOnode) then 209 write(6,'(/,a)') & 210 & 'siesta2wannier90: Reading the ' // trim(seedname) // '.nnkp file' 211 endif 212 call read_nnkp( seedname, latvec, reclatvec, numkpoints, & 213 kpointsfrac, nncount, nnlist, nnfolding, & 214 numproj, projections, numexcluded, excludedbands, & 215 w90_write_amn ) 216 217! Compute the vectors that connect each mesh k-point to its nearest neighbours 218 call chosing_b_vectors( kpointsfrac, nncount, nnlist, nnfolding, & 219 bvectorsfrac ) 220 221! Compute the number of bands for wannierization 222 call number_bands_wannier( numbandswan ) 223 numbands(ispin) = numbandswan(ispin) 224 225! Chose which bands are excluded from the wannierization procedure 226 call set_excluded_bands( ispin, numexcluded, excludedbands, numbands, & 227 & isexcluded, isincluded, numincbands ) 228 229! Compute the matrix elements of the plane wave, 230! for all the wave vectors that connect a given k-point to its nearest 231! neighbours 232 call compute_pw_matrix( nncount, bvectorsfrac ) 233 234! Compute the coefficients of the wavefunctions at the 235! k-points of the Wannier90 mesh 236 call diagonalizeHk( ispin ) 237 238! Compute the overlap between periodic parts of the wave functions 239 if( w90_write_mmn ) call mmn( ispin ) 240 241! Compute the overlaps between Bloch states onto trial localized orbitals 242 if( w90_write_amn ) call amn( ispin ) 243 244! Write the eigenvalues in a file, in the format required by Wannier90 245 if( IOnode .and. w90_write_eig ) call writeeig( ispin ) 246 247! Write the values of the Bloch states in a box 248 if( w90_write_unk ) call writeunk( ispin ) 249 250 enddo 251 252 if (IOnode) then 253 write(6,'(/,a)') & 254 & 'siesta2wannier90: All the information dumped in the corresponding files' 255 write(6,'(a)') & 256 & 'siesta2wannier90: End of the interface between Siesta and Wannier90' 257 endif 258 259 call timer("siesta2wannier90",2) 260 261end subroutine siesta2wannier90 262 263function getFileNameRoot(ispin,nspin,root) 264! 265! Constructs filenames for various cases of spin: 266! Nonpolarized, down and up by adding an extension to 267! "root". Behavior: 268! (a) unpolarized: no change 269! (b) spin up: _up suffix 270! (c) spin down: _dn suffix 271! 272 character(*),intent(in) :: root 273 integer,intent(in) :: ispin, nspin 274 character(len_trim(root)+3) :: getFileNameRoot 275 276 select case (nspin) 277! If only one spin component, do nothing 278 case (1) 279 getFileNameRoot = trim(root) 280! If two spin components, 281! append "_up" (for the first one) or "_dn" (for the second spin component) 282! to the seed name 283 case (2) 284 select case (ispin) 285 case (1) 286 getFileNameRoot = trim(root)//"_up" 287 case (2) 288 getFileNameRoot = trim(root)//"_dn" 289 end select 290! If more than two spin components, stop the program 291! and print and error message 292 case default 293 call die("getFileNameRoot: non-collinear spin not implemented yet") 294 end select 295end function getFileNameRoot 296 297endmodule m_siesta2wannier90 298