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! --- 8!================================================================== 9! 10program fatband 11! 12! 13! Computes the generalized projection of each eigenvector on a given 14! orbital set. 15! The syntax for orbital sets is the same as in mprop, so this 16! program can re-use as description file the same DOS-type file 17! as mprop. 18! 19! 20 21 use main_vars 22 use orbital_set, only: get_orbital_set 23 use io_hs, only: read_hs_file 24 use read_curves, only: read_curve_information, mask_to_arrays 25 26 implicit none 27 28 logical :: gamma_wfsx, got_qcos 29 integer :: ii1, ii2, ind, ind_red, no1, no2, n_int, nnz 30 real(dp) :: factor 31 32 real(dp), parameter :: tol_overlap = 1.0e-10_dp 33 34 logical, allocatable :: mask2(:) 35 integer, allocatable :: num_red(:), ptr(:), list_io2(:), list_ind(:) 36 real(dp), allocatable :: eig(:,:), fat(:,:) 37 38 39 integer :: nwfmx, nwfmin 40 integer :: min_band = 1 41 integer :: max_band = huge(1) 42 logical :: min_band_set = .false. 43 logical :: max_band_set = .false. 44 logical :: band_interval_set = .false. 45 real(dp) :: min_eigval 46 real(dp) :: max_eigval 47 real(dp) :: min_eigval_in_file 48 real(dp) :: min_eigval_in_band_set 49 real(dp) :: max_eigval_in_band_set 50 real(dp) :: minimum_spec_eigval = -huge(1.0_dp) 51 real(dp) :: maximum_spec_eigval = huge(1.0_dp) 52 53 integer :: ib, nbands, nspin_blocks 54 logical :: non_coll 55 56 ! 57 ! Process options 58 ! 59 n_opts = 0 60 do 61 call getopts('dhls:n:m:M:R:b:B:',opt_name,opt_arg,n_opts,iostat) 62 if (iostat /= 0) exit 63 select case(opt_name) 64 case ('d') 65 debug = .true. 66 case ('l') 67 energies_only = .true. 68 case ('R') 69 ref_line_given = .true. 70 ref_line = opt_arg 71 case ('b') 72 read(opt_arg,*) min_band 73 min_band_set = .true. 74 case ('B') 75 read(opt_arg,*) max_band 76 max_band_set = .true. 77 case ('h') 78 call manual() 79 case ('?',':') 80 write(0,*) "Invalid option: ", opt_arg(1:1) 81 write(0,*) "Usage: fat [ options ] DescriptionFile " 82 write(0,*) "Use -h option for manual" 83 STOP 84 end select 85 enddo 86 87 nargs = command_argument_count() 88 nlabels = nargs - n_opts + 1 89 if (nlabels /= 1) then 90 write(0,*) "Usage: fat [ options ] DescriptionFile " 91 write(0,*) "Use -h option for manual" 92 STOP 93 endif 94 95 call get_command_argument(n_opts,value=mflnm,status=iostat) 96 if (iostat /= 0) then 97 STOP "Cannot get .mpr file root" 98 endif 99 100 band_interval_set = (min_band_set .or. max_band_set) 101 102 !================================================== 103 104 ierr=0 105 106 ! Read type of job 107 108 open(mpr_u,file=trim(mflnm) // ".mpr", status='old') 109 read(mpr_u,*) sflnm 110 read(mpr_u,*) what 111 if (trim(what) /= "DOS") STOP "Fatbands needs DOS-style jobfile" 112 113 !================================================== 114 ! Read WFSX file 115 116 write(6,"(a)") "Reading wf file: " // trim(sflnm) // ".WFSX" 117 118 open(wfs_u,file=trim(sflnm)//'.WFSX',status='old',form='unformatted') 119 read(wfs_u) nkp, gamma_wfsx 120 allocate (wk(nkp), pk(3,nkp)) 121 122 read(wfs_u) nsp 123 non_coll = (nsp >= 4) 124 read(wfs_u) nao 125 read(wfs_u) !! Symbols, etc 126 if (debug) print *, "WFSX read: nkp, nsp, nnao: ", nkp, nsp, nao 127 128 nwfmx = -huge(1) 129 nwfmin = huge(1) 130 min_eigval = huge(1.0_dp) 131 max_eigval = -huge(1.0_dp) 132 min_eigval_in_band_set = huge(1.0_dp) 133 max_eigval_in_band_set = -huge(1.0_dp) 134 135 if (non_coll) then 136 nspin_blocks = 1 137 else 138 nspin_blocks = nsp 139 endif 140 141 do ik=1,nkp 142 do is=1,nspin_blocks 143 144 read(wfs_u) idummy, pk(1:3,ik), wk(ik) 145 if (idummy /= ik) stop "ik index mismatch in WFS file" 146 read(wfs_u) is0 147 read(wfs_u) number_of_wfns 148 nwfmx = max(nwfmx,number_of_wfns) 149 nwfmin = min(nwfmin,number_of_wfns) 150 151 do iw=1,number_of_wfns 152 read(wfs_u) iw0 153 read(wfs_u) eigval 154 min_eigval = min(min_eigval,eigval) 155 max_eigval = max(max_eigval,eigval) 156 ! 157 ! 158 if ((iw>=min_band).and.(iw<=max_band)) then 159 min_eigval_in_band_set = min(min_eigval_in_band_set,eigval) 160 max_eigval_in_band_set = max(max_eigval_in_band_set,eigval) 161 endif 162 read(wfs_u) 163 enddo 164 enddo 165 enddo 166 167 print "(a,2i5)", " Minimum/Maximum number of wfs per k-point: ", nwfmin, nwfmx 168 print "(a,2f12.4)", "Min_eigval, max_eigval on WFS file: ", & 169 min_eigval, max_eigval 170 print "(a,2f12.4)", "Min_eigval, max_eigval in band set : ", & 171 min_eigval_in_band_set, max_eigval_in_band_set 172 173 if (.not. band_interval_set) then 174 175 min_band = 1 ! Already set by default 176 max_band = nwfmin 177 178 else 179 180 if (min_band_set .and. (min_band < 1)) then 181 print "(a)", " ** Min_band implicitly reset to 1..." 182 min_band = 1 183 endif 184 if (min_band_set .and. (min_band > nwfmin)) then 185 print "(a,2i5)", " ** Min_band is too large for some k-points: (min_band, nwfmin):", min_band, nwfmin 186 STOP 187 endif 188 if (max_band_set .and. (max_band > nwfmin)) then 189 print "(a,2i5)", " ** Max_band is too large for some k-points: (max_band, nwfmin):", max_band, nwfmin 190 print "(a)", " ** Max_band will be effectively reset to its maximum allowed value" 191 max_band = nwfmin 192 endif 193 if (max_band_set .and. (max_band < min_band)) then 194 print "(a,2i5)", " ** Max_band is less than min_band: (max_band, min_band):", max_band, min_band 195 STOP 196 endif 197 198 min_eigval = min_eigval_in_band_set 199 max_eigval = max_eigval_in_band_set 200 endif 201 print "(a,3i4)", "Band set used: (min, max):", min_band, max_band 202 203 ! Read HSX file 204 ! Will pick up atoms, zval, and thus the nominal number of electrons, 205 ! but the total charge is read as qtot. 206 207 call read_hs_file(trim(sflnm)//".HSX") 208 209 if (energies_only) STOP 210 211 212 !==================== 213 214 ! * Orbital list 215 216 allocate(za(no_u), zc(no_u), zn(no_u), zl(no_u), zx(no_u), zz(no_u)) 217 nao = 0 218 do ia=1,na_u 219 it = isa(ia) 220 io = 0 221 do 222 io = io + 1 223 if (io > no(it)) exit 224 lorb = lquant(it,io) 225 do ko = 1, 2*lorb + 1 226 nao = nao + 1 227 za(nao)=ia 228 zc(nao)=it 229 zn(nao)=nquant(it,io) 230 zl(nao)=lorb 231 zx(nao)=ko 232 zz(nao)=zeta(it,io) 233 enddo 234 io = io + 2*lorb 235 enddo 236 enddo 237 if (nao /= no_u) STOP "nao /= no_u" 238 239 ! ================================== 240 241! 242! Process orbital sets 243! 244 allocate(orb_mask(no_u,2,ncbmx)) 245 allocate (koc(ncbmx,2,no_u)) 246 247 ! Give values to flags for reuse of the reading routine 248 dos = .true. 249 coop = .false. 250 call read_curve_information(.true.,.false., & 251 mpr_u,no_u,ncbmx,ncb,tit,orb_mask,dtc) 252 253 orb_mask(:,2,1:ncb) = .true. ! All orbitals considered 254 255 call mask_to_arrays(ncb,orb_mask(:,1,:),noc(:,1),koc(:,1,:)) 256 call mask_to_arrays(ncb,orb_mask(:,2,:),noc(:,2),koc(:,2,:)) 257 258!! 259 write(6,"('Writing files: ',a,'.stt ...')") trim(mflnm) 260 open(stt_u,file=trim(mflnm)//'.info') 261 write(stt_u,"(/'UNIT CELL ATOMS:')") 262 write(stt_u,"(3x,i4,2x,i3,2x,a20)") (i, isa(i), label(isa(i)), i=1,na_u) 263 write(stt_u,"(/'BASIS SET:')") 264 write(stt_u,"(5x,a20,3(3x,a1))") 'spec', 'n', 'l', 'z' 265 do it=1,nspecies 266 write(stt_u,"(5x,a20)") trim(label(it)) 267 io = 0 268 do 269 io = io + 1 270 if (io > no(it)) exit 271 write(stt_u,"(3(2x,i2))") nquant(it,io), lquant(it,io), zeta(it,io) 272 io = io + 2*lquant(it,io) 273 enddo 274 enddo 275 276 if ( nsp == 8 ) then 277 write(stt_u,"(/'SPIN (spin-orbit): ',i2)") nspin_blocks 278 else if ( nsp == 4 ) then 279 write(stt_u,"(/'SPIN (non-coll): ',i2)") nspin_blocks 280 else 281 write(stt_u,"(/'SPIN: ',i2)") nspin_blocks 282 endif 283 284 write(stt_u,"(/'AO LIST:')") 285 taux=repeat(' ',len(taux)) 286 do io=1,no_u 287 taux(1:30)=repeat(' ',30) 288 ik=1 289 if (zl(io).eq.0) then 290 taux(ik:ik)='s' 291 ik=0 292 elseif (zl(io).eq.1) then 293 if (zx(io).eq.1) taux(ik:ik+1)='py' 294 if (zx(io).eq.2) taux(ik:ik+1)='pz' 295 if (zx(io).eq.3) taux(ik:ik+1)='px' 296 ik=0 297 elseif (zl(io).eq.2) then 298 if (zx(io).eq.1) taux(ik:ik+2)='dxy' 299 if (zx(io).eq.2) taux(ik:ik+2)='dyz' 300 if (zx(io).eq.3) taux(ik:ik+2)='dz2' 301 if (zx(io).eq.4) taux(ik:ik+2)='dxz' 302 if (zx(io).eq.5) taux(ik:ik+5)='dx2-y2' 303 ik=0 304 elseif (zl(io).eq.3) then 305 taux(ik:ik)='f' 306 elseif (zl(io).eq.4) then 307 taux(ik:ik)='g' 308 elseif (zl(io).eq.5) then 309 taux(ik:ik)='h' 310 endif 311 write(stt_u,"(3x,i5,2x,i3,2x,a20)",advance='no') & 312 io, za(io), trim(label(zc(io))) 313 if (ik.eq.0) then 314 write(stt_u,"(3x,i2,a)") zn(io), trim(taux) 315 else 316 write(stt_u,"(3x,i2,a,i2.2)") zn(io), trim(taux), zx(io) 317 endif 318 enddo 319 320 write(stt_u,"(/'KPOINTS:',i7)") nkp 321 do ik=1,nkp 322 write(stt_u,"(3x,3f9.6)") pk(:,ik) 323 enddo 324 325 write(stt_u,"(/'FATBAND ORBITAL SETS:')") 326 do ic=1,ncb 327 write(stt_u,"(3x,a)") trim(tit(ic)) 328 write(stt_u,"(3x,'AO set I: ',/,15x,12i5)") (koc(ic,1,j),j=1,noc(ic,1)) 329 write(stt_u,"(3x,'Number of set II orbs:',i8)") noc(ic,2) 330 enddo 331 332 close(stt_u) 333 334 !================================== 335 336 if (ref_line_given) then 337 allocate(ref_mask(no_u)) 338 print *, "Orbital set spec: ", trim(ref_line) 339 call get_orbital_set(ref_line,ref_mask) 340 do io=1, no_u 341 if (ref_mask(io)) write(6,fmt="(i5)",advance="no") io 342 enddo 343 deallocate(ref_mask) 344 write(6,*) 345 STOP "bye from ref_line processing" 346 endif 347 348 349 !===================== 350 351 ! * Fatband weights 352 353 ! nspin has been read in iohs 354 nbands = max_band - min_band + 1 355 allocate(eig(nbands,nspin_blocks), fat(nbands,nspin_blocks)) 356 357 ! The first dimension is the number of real numbers per orbital 358 ! 1 for real wfs, 2 for complex, and four for the two spinor components 359 360 if (non_coll) then 361 allocate(wf_single(4,1:no_u)) 362 allocate(wf(4,1:no_u)) 363 else 364 if (gamma_wfsx) then 365 allocate(wf_single(1,1:no_u)) 366 allocate(wf(1,1:no_u)) 367 else 368 allocate(wf_single(2,1:no_u)) 369 allocate(wf(2,1:no_u)) 370 endif 371 endif 372 allocate (mask2(1:no_u)) 373 374 do ic=1,ncb 375 376 no1 = noc(ic,1) 377 no2 = noc(ic,2) 378 379 mask2(1:no_u) = .false. 380 do i2=1,no2 381 io2=koc(ic,2,i2) ! AO Set II 382 mask2(io2) = .true. 383 enddo 384 385 ! Create reduced pattern 386 ! First pass for checking dimensions 387 388 allocate (num_red(no1)) 389 do i1=1,no1 390 num_red(i1) = 0 391 io1=koc(ic,1,i1) ! AO Set I 392 do ii1 = 1,numh(io1) 393 ind = listhptr(io1)+ii1 394 ii2 = indxuo(listh(ind)) ! Equiv orb in unit cell 395 396 if ( .not. mask2(ii2)) Cycle ! Is not one of Set II 397 398 ! No distance restrictions for PDOS or FATBANDS, just an overlap check 399 400 if ( abs(Sover(ind)) < tol_overlap) CYCLE ! Not overlapping 401 402 num_red(i1) = num_red(i1) + 1 403 enddo 404 enddo 405 allocate (ptr(no1)) 406 ptr(1)=0 407 do i1=2,no1 408 ptr(i1)=ptr(i1-1)+num_red(i1-1) 409 enddo 410 nnz = sum(num_red(1:no1)) 411 412 write(*,"(a,3x,a,2x,a,i6,1x,i12)") 'Fatband coeffs set: ', trim(tit(ic)), & 413 'Base orbitals and interactions: ', & 414 no1, nnz 415 416 allocate (list_io2(nnz)) 417 allocate (list_ind(nnz)) 418 419 n_int = 0 420 do i1=1,no1 421 io1=koc(ic,1,i1) ! AO Set I 422 do ii1 = 1,numh(io1) 423 ind = listhptr(io1)+ii1 424 ii2 = indxuo(listh(ind)) 425 426 if ( .not. mask2(ii2)) Cycle 427 428 ! No distance restrictions for PDOS, just an overlap check 429 430 if ( abs(Sover(ind)) < tol_overlap) CYCLE ! Not overlapping 431 432 n_int = n_int + 1 433 list_io2(n_int) = ii2 434 list_ind(n_int) = ind 435 enddo 436 enddo 437 if (n_int .ne. nnz) then 438 print *, "n_int, nnz:", n_int, nnz 439 STOP "mismatch" 440 endif 441 442 443 !Stream over file, without using too much memory 444 445 rewind(wfs_u) 446 447 read(wfs_u) 448 read(wfs_u) 449 read(wfs_u) 450 read(wfs_u) 451 452 open(fat_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.EIGFAT') 453 write(fat_u,"(a,2i5)") "# " // trim(sflnm) // " min_band, max_band: ", min_band, max_band 454 write(fat_u,"(3i6)") nbands, nspin_blocks, nkp 455 456 do ik=1,nkp 457 458 write(fat_u,"(i4,3(1x,f10.5))") ik, pk(1:3,ik) 459 460 do is=1,nspin_blocks 461 462 ib = 0 463 fat(:,is) = 0.0_dp 464 465 read(wfs_u) 466 read(wfs_u) 467 read(wfs_u) number_of_wfns 468 do iw=1,number_of_wfns 469 read(wfs_u) 470 read(wfs_u) eigval 471 472 ! Use only the specified band set 473 if ( (iw<min_band) .or. (iw>max_band)) then 474 read(wfs_u) ! Still need to read this 475 CYCLE 476 endif 477 478 ib = ib + 1 479 eig(ib,is) = eigval ! This will be done for every curve... harmless 480 481 read(wfs_u) (wf_single(:,io), io=1,no_u) 482 ! Use a double precision form in what follows 483 wf(:,:) = real(wf_single(:,:), kind=dp) 484 485 486 do i1 = 1, no1 487 io1=koc(ic,1,i1) ! AO Set I 488 489 do i2 = 1,num_red(i1) 490 ind_red = ptr(i1)+i2 491 io2 = list_io2(ind_red) 492 ind = list_ind(ind_red) 493 494 ! (qcos, qsin) = C_1*conjg(C_2) 495 !AG: Corrected: (qcos, qsin) = conjg(C_1)*(C_2) 496 ! We might want to avoid recomputing this 497 if (non_coll) then 498 ! Take as weight the "complete spinor" product 499 qcos= wf(1,io1)*wf(1,io2) + & 500 wf(2,io1)*wf(2,io2) + & 501 wf(3,io1)*wf(3,io2) + & 502 wf(4,io1)*wf(4,io2) 503 qsin= wf(1,io1)*wf(2,io2) - & 504 wf(2,io1)*wf(1,io2) + & 505 wf(3,io1)*wf(4,io2) - & 506 wf(4,io1)*wf(3,io2) 507 else 508 if (gamma_wfsx) then 509 qcos = wf(1,io1)*wf(1,io2) 510 qsin = 0.0_dp 511 else 512 qcos= (wf(1,io1)*wf(1,io2) + & 513 wf(2,io1)*wf(2,io2)) 514 qsin= (wf(1,io1)*wf(2,io2) - & 515 wf(2,io1)*wf(1,io2)) 516 endif 517 endif 518 ! k*R_12 (r_2-r_1) 519 alfa=dot_product(pk(1:3,ik),xij(1:3,ind)) 520 521 ! Crb = Real(C_1*conjg(C_2)*exp(-i*alfa)) * S_12 522 !AG: This one better -- or Real(conjg(C_1)*C_2)*exp(+i*alfa)) * S_12 523 ! Common factor computed here 524 factor = (qcos*cos(alfa)-qsin*sin(alfa)) 525 526 fat(ib,is) = fat(ib,is) + Sover(ind)*factor 527 528 enddo ! i2 529 enddo ! i1 530 531 enddo ! iwf 532 write(fat_u,"(4(4x,f10.4,f9.5))") & 533 (eig(ib,is),fat(ib,is),ib=1,nbands) 534 enddo ! is 535 536 enddo ! ik 537 538 deallocate (num_red) 539 deallocate (ptr) 540 deallocate (list_io2) 541 deallocate (list_ind) 542 543 close(fat_u) 544 545 enddo ! ic 546 547 CONTAINS 548 549 subroutine manual() 550 551 write(6,"('* FAT(BANDS) PROGRAM')") 552 write(6,"(' Alberto Garcia, ICMAB-CSIC, 2012 ')") 553 write(6,*) 554 write(6,"(' FAT calculates eigenvector projections ')") 555 write(6,"(' using output files obtained with SIESTA. The atomic orbital (AO)')") 556 write(6,"(' sets are defined in an input file (MLabel.mpr).')") 557 write(6,"(' ')") 558 write(6,*) "Usage: fat [ options ] MPROP_FILE_BASENAME" 559 write(6,*) "Options:" 560 write(6,*) " -h: print manual " 561 write(6,*) " -d: debug " 562 write(6,*) " -l: print summary of energy information " 563 write(6,*) " " 564 write(6,*) " Selection of eigenstates to be used: " 565 write(6,*) " " 566 write(6,*) " -b Min_band : set minimum band index to be used " 567 write(6,*) " -B Max_band : set maximum band index to be used " 568 write(6,*) " " 569 write(6,*) 570 write(6,"('* .mpr FILE STRUCTURE')") 571 write(6,"(' SLabel # Name of the siesta output files')") 572 write(6,"(' DOS # DOS option of mprop is mandatory')") 573 write(6,"(' /- As many blocks as projections wanted ]')") 574 write(6,"(' | projection_name # DOS projection name')") 575 write(6,"(' \- Subset of AO (*) # Subset of orbitals included')") 576 write(6,"(' (*) See below how to define subsets of AO')") 577 write(6,"(' A final line with leading chars ---- can signal the end of the input')") 578 write(6,*) 579 write(6,"('* INPUT FILES')") 580 write(6,"(' [output files from SIESTA >= 2.4.1]')") 581 write(6,"(' SLabel.WFSX and SLabel.HSX (new format)')") 582 write(6,*) 583 write(6,"('* OUTPUT FORMAT')") 584 write(6,*) 585 write(6,*) " MLabel.CurveName.EIGFAT : File with eigenvalue and projection info " 586 write(6,"(' [An information file with extension .info will always be generated]')") 587 write(6,*) 588 write(6,"('* PROJECTION AND CURVES NAMES')") 589 write(6,"(' Alphanumerical string up to 30 char. with no spaces')") 590 write(6,"('* SUBSET OF AO USING ORDER NUMBERS')") 591 write(6,"(' List of integer numbers preceeded by a + symbol')") 592 write(6,"(' Each number refers to one AO in the final list of AO of SIESTA')") 593 write(6,"(' Example: + 23 65 78')") 594 write(6,"('* SUBSET OF AO USING ATOM_SHELL NOTATION')") 595 write(6,"(' List of atoms and shell groups of AO')") 596 write(6,"(' General notation: ATOM_SHELL')") 597 write(6,"(' > ATOM: Atomic symbol refers to all the atoms of that type')") 598 write(6,"(' Integer number refers to the N-th atom in unit cell')") 599 write(6,"(' > SHELL: Integer1+Letter+Integer2')") 600 write(6,"(' > Integer1 refers to the n quantum number')") 601 write(6,"(' > Letter refers to the l quantum number (s,p,d,f,g,h)')") 602 write(6,"(' > Integer2 refers to a single AO into the n-l shell')") 603 write(6,"(' Alternatively, alphanumerical strings can be used')") 604 write(6,"(' p-shells 1 y d-shells 1 xy 4 xz')") 605 write(6,"(' 2 z 2 yz 5 x2-y2')") 606 write(6,"(' 3 x 3 z2')") 607 write(6,"(' Particular cases:')") 608 write(6,"(' > Just ATOM is indicated: all the AO of the atom will be included')") 609 write(6,"(' > No value for Integer2: all the AO of the shell will be included')") 610 write(6,"(' Example: Ca_3p Al 4_4d3 5 O_2py')") 611 stop 612 613 end subroutine manual 614 615 616end program fatband 617 618 619