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_digest_nnkp 9 10 use precision, only: dp 11 12 implicit none 13 14 integer :: iun_nnkp ! Integer that points to the .nnkp file 15 ! used in read_nnkp and scan_file_to 16 17CONTAINS 18! subroutine read_nnkp 19! subroutine scan_file_to 20! function vectorproduct 21! subroutine chosing_b_vectors 22! subroutine set_excluded_bands 23! function getdelkmatgenhandle 24! subroutine number_bands_wannier 25 26 27subroutine read_nnkp( seedname, latvec, reclatvec, numkpoints, & 28 kpointsfrac, nncount, nnlist, nnfolding, & 29 numproj, projections, numexcluded, excludedbands, & 30 w90_write_amn ) 31!----------------------------------------------------------------------- 32! This subroutine reads the .nnkp file generated by the Wannier90 code 33! run in preprocessor mode, and populates some variables defined in the 34! siesta2wannier module. 35! The meaning of all these variables are properly described 36! in the file siesta2wannier90.F90 file 37!----------------------------------------------------------------------- 38 use sys, only: die ! Termination routine 39 use siesta_geom, only: ucell ! Unit cell lattice vectors 40 use units, only: Ang ! Conversion factor 41 use alloc, only: re_alloc ! Reallocation routines 42 use m_mpi_utils, only: broadcast ! Broadcasting routines 43 use files, only: label_length ! Number of characters allowed 44 ! in the name of a file 45 use parallel, only: IOnode ! Input/output node 46 use parallel, only: Node ! Working node 47 use trialorbitalclass ! Derived class to define the 48 ! trial wave functions 49 50 implicit none 51 52 character(label_length+3), intent(in) :: seedname ! Name of input file 53 integer, intent(out) :: numkpoints 54 integer, intent(out) :: nncount 55 integer, intent(out) :: numproj 56 integer, intent(out) :: numexcluded 57 integer, pointer :: nnlist(:,:) 58 integer, pointer :: nnfolding(:,:,:) 59 integer, pointer :: excludedbands(:) 60 real(dp), intent(out) :: latvec(3,3) 61 real(dp), intent(out) :: reclatvec(3,3) 62 real(dp), pointer :: kpointsfrac(:,:) 63 type(trialorbital), allocatable, intent(out) :: projections(:) 64 logical, intent(out) :: w90_write_amn 65 66 external :: io_assign, io_close 67 68! 69! Internal variables 70! 71 integer :: i, j, ik, iw! Loop counters 72 integer :: ikp, inn ! Loop counters 73 integer :: iikp ! Dummy variable 74 integer :: iu ! Integer number to refer to the nnkp file 75 logical :: have_nnkp ! Yes if the .nnkp is found 76 ! No if the .nnkp file is not found 77 real(dp) :: rcell(3,3) ! Reciprocal lattice vectors 78! Variables related to the projection functions. 79 real(dp) :: center_w(3) ! Projection function center, 80 ! as read from the nnkp file 81 ! (relative to the direct lattice 82 ! vectors) 83 ! They are transformed latter to 84 ! cartesian coordinates. 85 integer :: l_w ! Angular momentum of the projection func. 86 integer :: mr_w ! z-projection quantum number of the 87 ! projection function 88 integer :: r_w ! Radial quantum number of the 89 ! projection function 90 real(dp) :: zaxis(3) ! Defines the axis from which the polar 91 ! angle theta in spherical 92 ! polar coordinates is measured. 93 ! Default: (0.0 0.0 1.0) 94 real(dp) :: xaxis(3) ! Defines the axis from which the 95 ! azimuthal angle phi in spherical 96 ! polar coordinates is measured. 97 ! Must be orthogonal to z-axis. 98 real(dp) :: zovera ! z/a, diffusivity, spread. 99 ! Units in nnkp file: Ang^-1. 100 ! This is transformed latter to Bohr^-1 101 real(dp) :: xnorm ! Norm of the x-vector 102 real(dp) :: znorm ! Norm of the z-vector 103 real(dp) :: coseno ! Cosine between the x and z vectors 104 105 real(dp), parameter :: eps4 = 1.0e-4_dp 106 real(dp), parameter :: eps6 = 1.0e-6_dp 107 108 109! Initialize some variables 110 nullify( kpointsfrac ) 111 nullify( nnlist ) 112 nullify( nnfolding ) 113 nullify( excludedbands ) 114 115! Check the existence of the .nnkp file and open it if it exists. 116 if (IOnode) then 117 118 inquire(file=trim(seedname)//".nnkp",exist=have_nnkp) 119 if( .not. have_nnkp ) then 120 write(6,'(/,a)') & 121 & 'read_nnkp: Could not find the file '//trim(seedname)//'.nnkp' 122 call die() 123 endif 124 125 call io_assign( iu ) 126 iun_nnkp = iu 127 open( iun_nnkp, file=trim(seedname)//".nnkp",form='formatted') 128 129 endif ! IOnode 130 131 if (IOnode) then 132! check the information from *.nnkp with the nscf_save data 133 write(6,'(/,a)') & 134 & 'read_nnkp: Checking info from the ' // trim(seedname) // '.nnkp file' 135 136! 137! Real lattice vectors 138! 139 call scan_file_to('real_lattice') 140 write(6,'(a)')'read_nnkp: Reading data about real lattice' 141 do j = 1, 3 142 read(iun_nnkp,*) (latvec(i,j),i=1,3) 143! Rescale from Ang to Bohr, that are the units internally used in Siesta 144 do i = 1, 3 145 latvec(i,j) = latvec(i,j)*Ang 146 enddo 147 enddo 148! Compare the unit cell read from .nnkp with the unit cell used in Siesta 149 do j = 1, 3 150 do i = 1, 3 151 if(abs(latvec(i,j)-ucell(i,j))>eps4) then 152 write(6,*) 'read_nnkp: Something wrong with the real lattice! ' 153 write(6,*) ' latvec(i,j) =',latvec(i,j),' ucell(i,j)=',ucell(i,j) 154 call die() 155 endif 156 enddo 157 enddo 158 write(6,'(a)')'read_nnkp: - Real lattice is ok' 159 160! 161! Reciprocal lattice vectors 162! 163 call scan_file_to('recip_lattice') 164 write(6,'(a)')'read_nnkp: Reading data about reciprocal lattice' 165 do j = 1, 3 166 read(iun_nnkp,*) (reclatvec(i,j),i=1,3) 167! Rescale from Ang^-1 to Bohr^-1, that are the units internally 168! used in Siesta 169 do i = 1, 3 170 reclatvec(i,j) = reclatvec(i,j)/Ang 171 enddo 172 enddo 173! Compute the reciprocal lattice from the unit cell used in Siesta 174 call reclat(ucell, rcell, 1) 175! Compare the reciprocal lattice read from .nnkp with the one used inSiesta 176 do j = 1, 3 177 do i = 1, 3 178 if(abs(reclatvec(i,j)-rcell(i,j))>eps4) then 179 write(6,*)'read_nnkp: Something wrong with the reciprocal lattice!' 180 write(6,*)' reclatvec(i,j)=',reclatvec(i,j), & 181 & ' rcell(i,j)=',rcell(i,j) 182 call die() 183 endif 184 enddo 185 enddo 186 write(6,'(a)')'read_nnkp: - Reciprocal lattice is ok' 187 endif ! IOnode 188! Broadcast the lattice vectors in real space and the reciprocal lattice vectors 189 call broadcast( latvec ) 190 call broadcast( reclatvec ) 191 192! 193! k-points 194! 195 196 if (IOnode) then 197 call scan_file_to('kpoints') 198 write(6,'(a)')'read_nnkp: Reading data about k-points' 199 read(iun_nnkp,*) numkpoints 200 endif 201 202! Broadcast information regarding the number of k-points 203! and allocate in all nodes the corresponding arrays containing information 204! about the k-points. 205 call broadcast( numkpoints ) 206 call re_alloc( kpointsfrac, 1, 3, 1, numkpoints, & 207 & name='kpointsfrac', routine='read_nnkp') 208 209 if (IOnode) then 210 do ik = 1, numkpoints 211 read(iun_nnkp,*) kpointsfrac(1:3,ik) 212 enddo 213 endif 214 215! Broadcast the information regarding the k-points 216 call broadcast( kpointsfrac ) 217 218! 219! Projections 220! 221 numproj = 0 222 if (IOnode) then ! Read nnkp file on ionode only 223 write(6,'(a)')'read_nnkp: Reading data about projection centers ' 224 call scan_file_to('projections') 225 read(iun_nnkp,*) numproj 226 endif ! IOnode 227 call broadcast( numproj ) 228 229 if( numproj .ne. 0 ) then 230 if (IOnode) then ! read from ionode only 231 allocate(projections(numproj)) 232 do iw = 1, numproj 233 read(iun_nnkp,*) (center_w(i), i=1,3), & 234 & l_w, & 235 & mr_w, & 236 & r_w 237 read(iun_nnkp,*) (zaxis(i), i = 1, 3), & 238 & (xaxis(i), i = 1, 3), & 239 & zovera 240 241! Check that the xaxis and the yaxis are orthogonal to each other 242 xnorm = sqrt(xaxis(1)*xaxis(1) + xaxis(2)*xaxis(2) + xaxis(3)*xaxis(3)) 243 if (xnorm < eps6) call die('read_nnkp: |xaxis| < eps ') 244 znorm = sqrt(zaxis(1)*zaxis(1) + zaxis(2)*zaxis(2) + zaxis(3)*zaxis(3)) 245 if (znorm < eps6) call die('read_nnkp: |zaxis| < eps ') 246 247 coseno = (xaxis(1)*zaxis(1) + xaxis(2)*zaxis(2) + xaxis(3)*zaxis(3)) / & 248 & xnorm/znorm 249 if (abs(coseno) > eps6) & 250 call die('read_nnkp: xaxis and zaxis are not orthogonal') 251 if (zovera < eps6) & 252 call die('read_nnkp: zovera value must be positive') 253 254! Now convert "center" from ScaledByLatticeVectors to Bohrs 255 do i=1,3 256 projections(iw)%center(i) = 0.0_dp 257 do j=1,3 258 projections(iw)%center(i) = center_w(j)*latvec(i,j) + & 259 & projections(iw)%center(i) 260 enddo 261 enddo 262 263 projections(iw)%zaxis = zaxis 264 projections(iw)%xaxis = xaxis 265! Compute the y-axis 266 projections(iw)%yaxis = vectorproduct(zaxis,xaxis) 267 projections(iw)%yaxis = projections(iw)%yaxis / & 268 & sqrt(dot_product(projections(iw)%yaxis,projections(iw)%yaxis)) 269 270 projections(iw)%zovera = zovera*Ang ! To Bohr^-1 271 projections(iw)%r = r_w 272 projections(iw)%l = l_w 273 projections(iw)%mr = mr_w 274 275! Select the maximum angular momentum 276 select case(l_w) 277 case(0:3) 278 projections(iw)%lmax = l_w 279 case(-3:-1) 280 projections(iw)%lmax = 1 ! spN hybrids 281 case(-5:-4) 282 projections(iw)%lmax = 2 ! spdN hybrids 283 case default 284 call die("Invalid l in read_nkp") 285 end select 286 287! Further checks 288 if ( mr_w .lt. 1 ) & 289 & call die("Invalid mr in read_nnkp") 290 if ( r_w .gt. 3 .or. r_w .lt. 1 ) then 291 call die("Invalid r_w in read_nnkp") 292 elseif ( l_w .ge. 0 .and. mr_w .gt. 2*l_w+1 ) then 293 call die("Invalid mr_w in read_nnkp") 294 elseif ( l_w .lt. 0 .and. mr_w .gt. 1-l_w ) then 295 call die("Invalid mr_w in read_nnkp") 296 endif 297 298! Cut-off initialization (in Bohr) 299 projections(iw)%rcut = cutoffs(r_w)/projections(iw)%zovera 300 301 enddo 302 endif ! IOnode 303 304! Broadcast the information about the projection functions 305 call broadcast_projections( ) 306 307 endif ! numproj ne 0 308 309!! For debugging 310! if (IOnode) then ! read from ionode only 311! write(6,*) 'Projections:' 312! do iw = 1, numproj 313! write(6,'(3f12.6,3i3,f12.6)') & 314! & projections(iw)%center(1:3), & 315! & projections(iw)%l, & 316! & projections(iw)%mr, & 317! & projections(iw)%r, & 318! & projections(iw)%zovera 319! enddo 320! endif ! IOnode 321!! End debugging 322 323 if ( numproj .eq. 0 ) then 324 if( IOnode ) then 325 write(6,'(/,a)')'read_nnkp: WARNING: No projections specified in ' & 326 & //trim(seedname)//".nnkp" 327 write(6,'(a,/)')'read_nnkp: Calculation of Amn omitted' 328 w90_write_amn = .false. 329 endif 330 endif 331 332! 333! k-points neighbours 334! 335 if (IOnode) then ! read from ionode only 336 write(6,'(a)')'read_nnkp: Reading data about k-point neighbours ' 337 call scan_file_to('nnkpts') 338 read(iun_nnkp,*) nncount 339 endif 340 341! Broadcast information regarding the number of k-points neighbours 342! and allocate in all nodes the corresponding arrays containing information 343! about k-point neighbours 344 call broadcast( nncount ) 345 call re_alloc( nnlist, 1, numkpoints, 1, nncount, & 346 & name = "nnlist", routine = "read_nnkp" ) 347 call re_alloc( nnfolding, 1, 3, 1, numkpoints, 1, nncount, & 348 & name = "nnfolding", routine = "read_nnkp" ) 349 350 if (IOnode) then 351 do ikp = 1, numkpoints 352 do inn = 1, nncount 353 read(iun_nnkp,*) iikp, & 354 & nnlist(ikp,inn), & 355 & ( nnfolding(j,ikp,inn), j = 1, 3 ) 356 enddo 357 enddo 358 endif 359 360! Broadcast information regarding the list of k-points neighbours 361 call broadcast( nnlist ) 362 call broadcast( nnfolding ) 363 364! 365! excluded bands 366! 367 if (IOnode) then ! read from ionode only 368 write(6,'(a)')'read_nnkp: Reading data about excluded bands ' 369 call scan_file_to('exclude_bands') 370 read (iun_nnkp,*) numexcluded 371 endif 372 373! Broadcast information regarding the number of excluded bands 374! and allocate in all nodes the corresponding arrays containing 375! information about excluded bands 376 call broadcast( numexcluded ) 377 if( numexcluded .gt. 0 ) then 378 call re_alloc( excludedbands, 1, numexcluded, & 379 & name = 'excludedbands', routine = "read_nnkp" ) 380 endif 381 382 if (IOnode) then 383 if( numexcluded .gt. 0 ) then 384 do i = 1, numexcluded 385 read(iun_nnkp,*) excludedbands(i) 386 enddo 387 endif 388 endif 389 390!! For debugging 391! if (IOnode) then 392! do i = 1, numexcluded 393! write(6,*) excludedbands(i) 394! enddo 395! endif 396!! End debugging 397 398! Broadcast information regarding excluded bands 399 if( numexcluded .gt. 0 ) call broadcast( excludedbands ) 400 401 if (IOnode) call io_close(iun_nnkp) ! ionode only 402 403 return 404 405end subroutine read_nnkp 406 407! 408!----------------------------------------------------------------------- 409subroutine scan_file_to (keyword) 410! 411! This subroutine is used to digest the .nnkp file. 412! It is borrowed from the QuantumEspresso suite. 413! 414! 415 implicit none 416 character(len=*) :: keyword 417 character(len=80) :: line1, line2 418! 419! by uncommenting the following line the file scan restarts every time 420! from the beginning thus making the reading independent on the order 421! of data-blocks 422! rewind (iun_nnkp) 423! 42410 continue 425 read(iun_nnkp,*,end=20) line1, line2 426 if(line1/='begin') goto 10 427 if(line2/=keyword) goto 10 428 return 42920 write (6,*) keyword," data-block missing " 430 call die("--") 431end subroutine scan_file_to 432!----------------------------------------------------------------------- 433 434! 435!----------------------------------------------------------------------- 436function vectorproduct(a,b) 437! 438! This function computes the vector product of two vectors a and b 439! 440 real(dp),intent(in),dimension(3) :: a,b 441 real(dp),dimension(3) :: vectorproduct 442 vectorproduct(1) = a(2) * b(3) - a(3) * b(2) 443 vectorproduct(2) = a(3) * b(1) - a(1) * b(3) 444 vectorproduct(3) = a(1) * b(2) - a(2) * b(1) 445end function vectorproduct 446!----------------------------------------------------------------------- 447 448! 449!----------------------------------------------------------------------- 450subroutine chosing_b_vectors( kpointsfrac, nncount, nnlist, nnfolding, & 451 bvectorsfrac ) 452! 453! In this subroutine we compute the b vectors that connect each 454! mesh k-point to its nearest neighbours. 455! 456! ------------------------------ INPUT ----------------------------------------- 457! real(dp) kpointsfrac(:,:) : List of k points where the overlap matrix 458! of the periodic part of the wave function 459! with its neighbour will be computed. 460! Given in reduced units. 461! integer nncount : The number of nearest neighbours belonging to 462! each k-point of the Monkhorst-Pack mesh 463! integer nnlist(:,:) : nnlist(ikp,inn) is the index of the 464! inn-neighbour of ikp-point 465! in the Monkhorst-Pack grid folded to the 466! first Brillouin zone 467! integer nnfolding(:,:,:) : nnFolding(i,ikp,inn) is the i-component 468! of the reciprocal lattice vector 469! (in reduced units) that brings 470! the inn-neighbour specified in nnlist 471! (which is in the first BZ) to the 472! actual \vec{k} + \vec{b} that we need. 473! In reciprocal lattice units. 474! ------------------------------ OUTPUT ---------------------------------------- 475! real(dp) bvectorsfrac(:,:) : Vectors that connect each mesh k-point 476! to its nearest neighbours. 477! ------------------------------ BEHAVIOUR ------------------------------------- 478! Set some kpointfrac, in this case the first in the list, and calculate 479! b's by subtracting it from it's neighbors. 480! They will be used as a handle to the "delkmat" matrix element of 481! a plane wave with wave vector b. 482! ------------------------------ UNITS ----------------------------------------- 483! bvectorsfrac are given in reduced units. 484! ------------------------------------------------------------------------------ 485 486! Used module procedures 487 use alloc, only: re_alloc ! Re-allocation routine 488 489 implicit none 490 491 integer, intent(in) :: nncount 492 integer, pointer :: nnlist(:,:) 493 integer, pointer :: nnfolding(:,:,:) 494 real(dp), pointer :: kpointsfrac(:,:) 495 real(dp), pointer :: bvectorsfrac(:,:) 496 497! 498! Internal variables 499! 500 integer :: i, inn 501 502 nullify( bvectorsfrac ) 503 call re_alloc( bvectorsfrac, 1, 3, 1, nncount, & 504 name="bvectorsfrac", routine = "chosing_b_vectors") 505! Set some kvector, in this case the first in the list and calculate 506! b's by subtracting it from it's neighbors 507! They will be used as a handle to the "delkmat" matrix element of 508! a plane wave with wave vector b. 509! 510! Loop on the three spatial directions 511 do i = 1, 3 512! Loop on neighbour k-points 513 do inn = 1, nncount 514 bvectorsfrac(i,inn) = & 515 kpointsfrac(i,nnlist(1,inn)) + nnfolding(i,1,inn) - kpointsfrac(i,1) 516 enddo 517 enddo 518 519end subroutine chosing_b_vectors 520 521 522! 523!----------------------------------------------------------------------- 524integer function getdelkmatgenhandle( vec, nncount, bvectorsfrac ) 525! 526! This function accepts a k vector connecting two neighboring 527! k-points and returns the position of vec in the array bVectorsFrac. 528! It is used to assign matrix elements of a planewave to a k-vector. 529! Vec must be in fractional coordinates. 530 531 implicit none 532 real(dp), intent(in) :: vec(3) 533 integer, intent(in) :: nncount ! Number of nearest k-point 534 ! neighbours 535 real(dp), intent(in) :: bvectorsfrac(3,nncount) ! Vectors that connect each 536 ! mesh k-point to its 537 ! nearest neighbours. 538 539! Internal variables 540 integer :: inn 541 real(dp) :: normSq,minNormSq 542 543 minNormSq = 1.0_dp 544 do inn = 1, nncount 545 normSq = (bvectorsfrac(1,inn) - vec(1))**2 + & 546 & (bvectorsfrac(2,inn) - vec(2))**2 + & 547 & (bvectorsfrac(3,inn) - vec(3))**2 548 if ( normsq .lt. minnormsq) then 549 getdelkmatgenhandle = inn 550 minnormsq = normsq 551 endif 552 enddo 553end function getdelkmatgenhandle 554 555! 556!----------------------------------------------------------------------- 557subroutine number_bands_wannier( numbandswan ) 558 559 use siesta_options, only: hasnobup ! Is the number of bands with spin up 560 ! for wannierization defined 561 use siesta_options, only: hasnobdown ! Is the number of bands with spin down 562 ! for wannierization defined 563 use siesta_options, only: hasnob ! Is the number of bands 564 ! for wannierization defined 565 ! (for non spin-polarized calculations) 566 use siesta_options, only: nobup ! Number of bands with spin up 567 ! for wannierization 568 use siesta_options, only: nobdown ! Number of bands with spin down 569 ! for wannierization 570 use siesta_options, only: nob ! Number of bands for wannierization 571 ! (for non spin-polarized calculations) 572 use m_spin, only: nspin ! Number of spin components 573 use sparse_matrices, only: maxnh ! Maximum number of orbitals 574 ! interacting 575 use sparse_matrices, only: numh ! Number of nonzero elements of each 576 ! row of the hamiltonian matrix 577 use sparse_matrices, only: listhptr ! Pointer to start of each row of the 578 ! hamiltonian matrix 579 use sparse_matrices, only: Dscf ! Density matrix in sparse form 580 use sparse_matrices, only: S ! Overlap matrix in sparse form 581 use atomlist, only: no_l ! Number of orbitals in local node 582 ! NOTE: When running in parallel, 583 ! this is core dependent 584 ! Sum_{cores} no_l = no_u 585 use atomlist, only: no_u ! Number of orbitals in unit cell 586 587 use m_noccbands, only: number_of_occupied_bands 588 use m_noccbands, only: noccupied ! Number of occupied bands for a given 589 ! spin direction 590 use parallel, only: IOnode 591 592#ifdef MPI 593 use m_mpi_utils, only: globalize_sum 594#endif 595 596 implicit none 597 598! integer numbandswan(2) ! Number of bands for wannierization 599 integer, intent(out), dimension(2) :: numbandswan 600 601! Internal variables 602 integer :: ispin ! Counter for the loops on the spin variables 603 integer :: io ! Counter for loops on the atomic orbitals 604 integer :: j ! Counter for loops on neighbours 605 integer :: ind ! Index of the neighbour orbital 606 607 real(dp):: qspin(2) ! Total population of spin up and down 608 609#ifdef MPI 610 real(dp):: qtmp(2) ! Temporary for globalized spin charge 611#endif 612 613! Initialize variables 614 numbandswan(:) = 0 615 616! Calculate the default number of bands considered for wannierization 617! numbandswan = number of occupied bands 618! This makes sense in molecules and insulators. 619 do ispin = 1, nspin 620 qspin(ispin) = 0.0_dp 621 do io = 1, no_l 622 do j = 1, numh(io) 623 ind = listhptr(io) + j 624 qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind) 625 enddo 626 enddo 627 enddo 628 629#ifdef MPI 630! Global reduction of spin components 631 call globalize_sum(qspin(1:nspin),qtmp(1:nspin)) 632 qspin(1:nspin) = qtmp(1:nspin) 633#endif 634 635!! For debugging 636! if( IOnode ) then 637! write(6,'(/,a,2f12.5)') & 638! & 'number_bands_wannier: qspin', qspin 639! endif 640!! End debugging 641 642 call number_of_occupied_bands( qspin ) 643 644! Setting up the default number of bands for wannierization to the 645! number of occupied bands 646 do ispin = 1, nspin 647 numbandswan(ispin) = noccupied(ispin) 648 enddo 649 650! Check the consistency with the required bands specified in the input file 651 if ( nspin .eq. 1 ) then ! Unpolarized calculation 652! If a number of bands is explicitly included in the input, 653! then use it overwriting the default value 654 if ( hasnob ) then 655 numbandswan(1) = nob 656! Since we are in a non-spin polarized calculation 657! ignore any possible value of the number of the number of bands for 658! wannierization introduced in the input 659 if ( hasnobup ) then 660 if( IOnode ) & 661 & write(6,'(a)') & 662 & 'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBandsUp' 663 endif 664 if ( hasnobdown ) then 665 if( IOnode ) & 666 & write(6,'(a)') & 667 & 'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBandsDown' 668 endif 669 670 elseif ( hasnobup .and. hasnobdown ) then 671 numbandswan(1) = nobup 672 numbandswan(2) = nobdown 673 if( numbandswan(1) .ne. numbandswan(2) ) then 674 if( IOnode ) & 675 & write(6,'(a)') & 676 & 'number_bands_wannier: inconsistent output band specification' 677 call die("Confusion in fdf") 678 endif 679 elseif ( hasnobup ) then 680 numbandswan(1) = nobup 681 elseif ( hasnobdown ) then 682 numbandswan(1) = nobdown 683 else 684 numbandswan(1) = noccupied(1) 685 if( IOnode ) & 686 & write(6,'(/,a)') & 687 & 'number_bands_wannier: Automatic number of bands for wannierization' 688 endif 689 if ( numbandswan(1) .gt. no_u ) & 690 & call die( 'More bands required than available number of states. & 691 & Reduce the number of bands.') 692 else ! Spin polarized calculation 693 if ( hasnob ) then 694 if ( .not. ( hasnobup .or. hasnobdown )) then 695 numbandswan(1) = nob 696 numbandswan(2) = nob 697 elseif ( .not. ( hasnobup .and. hasnobdown )) then 698 if( IOnode ) & 699 & write(6,'(a)') & 700 & 'number_bands_wannier: inconsistent output band specification' 701 call die('Confusion in fdf') 702 else 703 if( IOnode ) & 704 & write(6,'(a)') & 705 & 'number_bands_wannier: ignoring Siesta2Wannier90.NumberOfBands' 706 numbandswan(1) = nobup 707 numbandswan(2) = nobdown 708 endif 709 elseif ( hasnobup .and. hasnobdown ) then 710 numbandswan(1) = nobup 711 numbandswan(2) = nobdown 712 elseif ( hasnobup .or. hasnobdown ) then 713 if( IOnode ) & 714 & write(6,'(a)') & 715 & 'number_bands_wannier: inconsistent output band specification' 716 call die('Confusion in fdf') 717 else 718 if( IOnode ) & 719 & write(6,'(/,a)') & 720 & 'number_bands_wannier: Automatic number of bands for wannierization' 721 endif 722 if ( numbandswan(1) .gt. no_u .or. numbandswan(2) .gt. no_u ) & 723 & call die( 'More bands required than available number of states. & 724 & Reduce the number of bands.') 725 endif 726 727 if( IOnode ) then 728 write(6,'(a)') & 729 & 'number_bands_wannier: Number of bands for wannierization ' 730 write(6,'(a,2i5)') & 731 & 'number_bands_wannier: before excluding bands = ', & 732 & numbandswan 733 endif 734 735end subroutine number_bands_wannier 736 737! 738!----------------------------------------------------------------------- 739subroutine set_excluded_bands( ispin, numexcluded, excludedbands, numbands, & 740 & isexcluded, isincluded, numincbands ) 741! 742! In this subroutine we set up the bands that are excluded from wannierization 743! 744! Used module procedures 745 use alloc, only: re_alloc ! Re-allocation routine 746 use alloc, only: de_alloc ! De-allocation routine 747 748! Used module parameters 749 use atomlist, only: no_u ! Number of orbitals in unit cell 750 use parallel, only: IOnode ! Input/output node 751 752!! For debugging 753 use parallel, only: Node ! Working node 754 use sys, only: die ! Termination routine 755!! End debugging 756 757 implicit none 758 759 integer, intent(in) :: ispin 760 integer, intent(in) :: numexcluded 761 ! Number of bands to exclude from the 762 ! from the calculation of the overlap 763 ! and projection matrices. 764 ! This variable is read from the .nnkp 765 ! file 766 integer, pointer :: excludedbands(:) 767 ! Bands to be excluded. 768 ! This variable is read from the .nnkp 769 ! file 770 integer, intent(in) :: numbands(2) 771 logical, pointer :: isexcluded(:) ! List of bands excluded for 772 ! Wannierization 773 integer, pointer :: isincluded(:) ! List of bands included for 774 ! Wannierization 775 integer, intent(out) :: numincbands(2) 776 777! 778! Internal variables 779! 780 integer :: iterex ! Counter for the loop on excluded bands 781 integer :: iband ! Index of the band that is excluded 782 783 nullify( isexcluded ) 784 call re_alloc( isexcluded, 1, no_u, name='isexcluded', & 785 & routine='set_excluded_bands' ) 786 787! By default, the occupied bands for the corresponding spin, numbands(ispin), 788! are not excluded from the calculation... 789 isexcluded(1:numbands(ispin)) = .false. 790 791! ... while the unoccupied bands (from numbands(ispin+1) to the total number 792! of bands (no_u)) are excluded 793 isexcluded(numbands(ispin)+1:no_u) = .true. 794 795 do iterex = 1, numexcluded 796 iband = excludedbands( iterex ) 797 if ( iband .lt. 1 ) then 798 call die("Erroneous list of excluded bands.") 799 elseif( iband .gt. numbands(ispin) ) then 800 call die("Erroneous list of excluded bands or numbands too small") 801 endif 802! Exclude the corresponding band from the computation 803 isexcluded( iband ) = .true. 804 enddo 805 806 call de_alloc( excludedbands, name='excludedbands', & 807 & routine='set_excluded_bands') 808 809! The number of actually "included" bands is the computation of 810! the overlap and projection matrices is the total number of 811! occupied bands minus the excluded bands 812! This number will be used to dimension Mmnkb and Amnmat matrices, 813! and also corrsponds to the number 814! of wave functions in UNK00001.1 815 numincbands(ispin) = numbands(ispin) - numexcluded 816 817! Determine the bands explicitly included for wannierization 818 nullify( isincluded ) 819 call re_alloc( isincluded, 1, numincbands(ispin), name='isincluded', & 820 & routine='set_excluded_bands' ) 821 822 iband = 1 823 do iterex = 1, numbands(ispin) 824 if ( isexcluded( iterex ) .eqv. .false. ) then 825 isincluded( iband ) = iterex 826 iband = iband + 1 827 endif 828 enddo 829 830 if( IOnode ) then 831 write(6,'(/,a,2i5)') 'Number of bands for wannierization ' // & 832 'after excluding bands:', numincbands 833 write(6,'(a)') 'Bands to be wannierized: ' 834 write(6,'(16i5)') (isincluded( iband ), iband=1, numincbands(ispin)) 835 endif 836 837end subroutine set_excluded_bands 838 839endmodule m_digest_nnkp 840