1! 2! Copyright (C) 2001-2012 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8PROGRAM plotband 9 10 ! reads data files produced by "bands.x", produces 11 ! * data file ready for plotting with gnuplot, xmgr or the like 12 ! * a postscript file that can be directly printed 13 ! for projected band structure, produces 14 ! * a gnuplot script that can be used to generate 15 ! postscript and/or pdf file 16 ! Important notice: 17 ! - k-points processed by bands.x should be along a continuous path 18 ! - no two consecutive k-points should be equal (i.e.: a k-point, 19 ! e.g. 0,0,0, can appear more than once but not in sequence) 20 ! If these rules are violated, unpredictable results may follow 21 22 !!! 23 USE parser 24 USE kinds, ONLY : DP 25 !!! 26 IMPLICIT NONE 27 INTEGER, PARAMETER :: stdout=6 28 real, ALLOCATABLE :: e(:,:), k(:,:), e_in(:), kx(:) 29 real :: k1(3), k2(3), ps 30 real, ALLOCATABLE :: e_rap(:,:), k_rap(:,:) 31 INTEGER, ALLOCATABLE :: nbnd_rapk(:), rap(:,:) 32 INTEGER, ALLOCATABLE :: npoints(:) 33 CHARACTER(len=1024) aux 34 INTEGER :: nks = 0, nbnd = 0, ios, nlines, n,i,j,ni,nf,nl, firstk, lastk 35 INTEGER :: nks_rap = 0, nbnd_rap = 0 36 LOGICAL, ALLOCATABLE :: high_symmetry(:), is_in_range(:), is_in_range_rap(:) 37 CHARACTER(len=256) :: filename, filename1, filenamegnu 38 NAMELIST /plot/ nks, nbnd 39 NAMELIST /plot_rap/ nks_rap, nbnd_rap 40 INTEGER :: n_interp 41 real, ALLOCATABLE :: k_interp(:), e_interp(:), coef_interp(:,:) 42 43 real :: emin = 1.e10, emax =-1.e10, etic, eref, deltaE, Ef 44 45 INTEGER, PARAMETER :: max_lines=99 46 real :: mine, dxmod, dxmod_save 47 INTEGER :: point(max_lines+1), nrap(max_lines) 48 INTEGER :: ilines, irap, ibnd, ipoint, jnow 49 50 real, PARAMETER :: cm=28.453, xdim=15.0*cm, ydim=10.0*cm, & 51 x0=2.0*cm, y0=2.0*cm, eps=1.e-4 52 53 LOGICAL :: exist_rap 54 LOGICAL, ALLOCATABLE :: todo(:,:) 55 CHARACTER(len=6), EXTERNAL :: int_to_char 56 !!! 57 LOGICAL :: exist_proj 58 CHARACTER(len=256) :: filename2, field 59 CHARACTER(len=:), ALLOCATABLE :: line 60 INTEGER :: nat, ntyp, natomwfc, nprojwfc, nwfc, idum 61 INTEGER, ALLOCATABLE :: atwfclst(:) 62 REAL(DP) :: proj, fdum 63 REAL(DP), ALLOCATABLE :: sumproj(:,:), p_rap(:,:) 64 !!! 65 66 67 CALL get_file ( filename ) 68 OPEN(unit=1,file=filename,form='formatted') 69 READ (1, plot, iostat=ios) 70 ! 71 IF (nks <= 0 .or. nbnd <= 0 .or. ios /= 0) THEN 72 STOP 'Error reading file header' 73 ELSE 74 WRITE(*,'("Reading ",i4," bands at ",i6," k-points")') nbnd, nks 75 ENDIF 76 77 filename1=trim(filename)//".rap" 78 !!! replace with inquire statement? 79 exist_rap=.true. 80 OPEN(unit=21,file=filename1,form='formatted',status='old',err=100,iostat=ios) 81 82100 IF (ios /= 0) THEN 83 exist_rap=.false. 84 ENDIF 85 !!! 86 IF (exist_rap) THEN 87 READ (21, plot_rap, iostat=ios) 88 IF (nks_rap/=nks.or.nbnd_rap/=nbnd.or.ios/=0) THEN 89 WRITE(6,'("file with representations not compatible with bands")') 90 exist_rap=.false. 91 ENDIF 92 ENDIF 93 ! 94 ALLOCATE (e(nbnd,nks)) 95 ALLOCATE (k(3,nks), e_in(nks), kx(nks), npoints(nks), high_symmetry(nks)) 96 ALLOCATE (is_in_range(nbnd)) 97 98 IF (exist_rap) THEN 99 ALLOCATE(nbnd_rapk(nks)) 100 ALLOCATE(e_rap(nbnd,nks)) 101 ALLOCATE(rap(nbnd,nks)) 102 ALLOCATE(k_rap(3,nks)) 103 ALLOCATE(todo(nbnd,2)) 104 ALLOCATE (is_in_range_rap(nbnd)) 105 ENDIF 106 !!! 107 filename2=trim(filename)//".proj" 108 exist_proj=.false. 109 INQUIRE(file=filename2,exist=exist_proj) 110 IF (exist_proj) THEN 111 OPEN(UNIT=22, FILE=filename2, FORM='formatted', STATUS='old', IOSTAT=ios) 112 IF (ios/=0) STOP 'Error opening projection file ' 113 READ (22, *, ERR=22, IOSTAT=ios) ! empty line 114 READ (22, '(8i8)', ERR=22, IOSTAT=ios) idum, idum, idum, idum, idum, & 115 idum, nat, ntyp ! FFT grid (ignored), nat, ntyp 116 READ (22, '(i6,6f12.8)', ERR=22, IOSTAT=ios) idum, fdum, fdum, fdum, fdum, & 117 fdum, fdum ! ibrav, celldm(1-6) 118 IF (idum == 0) THEN ! read and discard the three cell vectors 119 READ (22, *, ERR=22, IOSTAT=ios) 120 READ (22, *, ERR=22, IOSTAT=ios) 121 READ (22, *, ERR=22, IOSTAT=ios) 122 ENDIF 123 DO i=1,1+nat+ntyp 124 READ(22, *, ERR=22, IOSTAT=ios) ! discard atomic positions 125 ENDDO 126 READ (22, '(3i8)',ERR=22, IOSTAT=ios) natomwfc, nks_rap, nbnd_rap 127 READ (22, *, ERR=22, IOSTAT=ios) ! discard another line 128 129 IF (nks_rap/=nks.or.nbnd_rap/=nbnd.or.ios/=0) THEN 130 WRITE(6,'("file with projections not compatible with bands")') 131 exist_proj=.false. 132 ELSE 133 ALLOCATE( sumproj(nbnd,nks) ) 134 ALLOCATE( p_rap(nbnd,nks) ) 135 ENDIF 136 ENDIF 137 !!! 138 139 140 high_symmetry=.false. 141 142 DO n=1,nks 143 READ(1,*,end=20,err=20) ( k(i,n), i=1,3 ) 144 READ(1,*,end=20,err=20) (e(i,n),i=1,nbnd) 145 IF (exist_rap) THEN 146 READ(21,*,end=20,err=20) (k_rap(i,n),i=1,3), high_symmetry(n) 147 READ(21,*,end=20,err=20) (rap(i,n),i=1,nbnd) 148 IF (abs(k(1,n)-k_rap(1,n))+abs(k(2,n)-k_rap(2,n))+ & 149 abs(k(3,n)-k_rap(3,n)) > eps ) THEN 150 WRITE(stdout,'("Incompatible k points in rap file")') 151 DEALLOCATE(nbnd_rapk) 152 DEALLOCATE(e_rap) 153 DEALLOCATE(rap) 154 DEALLOCATE(k_rap) 155 DEALLOCATE(todo) 156 DEALLOCATE(is_in_range_rap) 157 CLOSE(unit=21) 158 exist_rap=.false. 159 ENDIF 160 ENDIF 161 ENDDO 162 CLOSE(unit=1) 163 IF (exist_rap) CLOSE(unit=21) 164 165 !!! 166 ! First read the list of atomic wavefunctions from the input, 167 ! then read the entire projection file and sum up the projections 168 ! onto the selected wavefunctions. 169 !!! 170 IF (exist_proj) THEN 171 WRITE(*,'("List of atomic wavefunctions: ")', advance="NO") 172 !read input string of arbitrary length 173 CALL readline(5,line) 174 CALL field_count( nprojwfc, line ) 175 ALLOCATE(atwfclst(nprojwfc)) 176 atwfclst(:) = -1 177 DO nwfc = 1,nprojwfc 178 CALL get_field(nwfc, field, line) 179 READ(field,*) atwfclst(nwfc) 180 ENDDO 181 182 sumproj(:,:) = 0.D0 183 DO nwfc = 1, natomwfc 184 READ(22, *, ERR=23, IOSTAT=ios) 185 DO n=1,nks 186 DO ibnd=1,nbnd 187 READ(22, '(2i8,f20.10)', ERR=23, IOSTAT=ios) idum,idum,proj 188 IF ( any( atwfclst(:) == nwfc ) ) sumproj(ibnd,n) = sumproj(ibnd,n) + proj 189 ENDDO 190 ENDDO 191 ENDDO 192 DEALLOCATE(atwfclst) 193 CLOSE(22) 194 ENDIF 195 !!! 196 197! 198! Now find the high symmetry points in addition to those already identified 199! in the representation file 200! 201 DO n=1,nks 202 IF (n==1 .or. n==nks) THEN 203 high_symmetry(n) = .true. 204 ELSE 205 k1(:) = k(:,n) - k(:,n-1) 206 k2(:) = k(:,n+1) - k(:,n) 207 ps = ( k1(1)*k2(1) + k1(2)*k2(2) + k1(3)*k2(3) ) / & 208 sqrt( k1(1)*k1(1) + k1(2)*k1(2) + k1(3)*k1(3) ) / & 209 sqrt( k2(1)*k2(1) + k2(2)*k2(2) + k2(3)*k2(3) ) 210 high_symmetry(n) = (abs(ps-1.d0) >1.0d-4).or.high_symmetry(n) 211! 212! The gamma point is a high symmetry point 213! 214 IF (k(1,n)**2+k(2,n)**2+k(3,n)**2 < 1.0d-9) high_symmetry(n)=.true. 215! 216! save the typical length of dk 217! 218 IF (n==2) dxmod_save = sqrt( k1(1)**2 + k1(2)**2 + k1(3)**2) 219 220 ENDIF 221 ENDDO 222 223 kx(1) = 0.d0 224 DO n=2,nks 225 dxmod=sqrt ( (k(1,n)-k(1,n-1))**2 + & 226 (k(2,n)-k(2,n-1))**2 + & 227 (k(3,n)-k(3,n-1))**2 ) 228 IF (dxmod > 5*dxmod_save) THEN 229! 230! A big jump in dxmod is a sign that the point k(:,n) and k(:,n-1) 231! are quite distant and belong to two different lines. We put them on 232! the same point in the graph 233! 234 kx(n)=kx(n-1) 235 ELSEIF (dxmod > 1.d-5) THEN 236! 237! This is the usual case. The two points k(:,n) and k(:,n-1) are in the 238! same path. 239! 240 kx(n) = kx(n-1) + dxmod 241 dxmod_save = dxmod 242 ELSE 243! 244! This is the case in which dxmod is almost zero. The two points coincide 245! in the graph, but we do not save dxmod. 246! 247 kx(n) = kx(n-1) + dxmod 248 249 ENDIF 250 ENDDO 251 252 DO n=1,nks 253 DO i=1,nbnd 254 emin = min(emin, e(i,n)) 255 emax = max(emax, e(i,n)) 256 ENDDO 257 ENDDO 258 WRITE(*,'("Range:",2f10.4,"eV Emin, Emax, [firstk, lastk] > ")', advance="NO") emin, emax 259 READ(5,'(a1024)') aux 260 READ(aux,*,iostat=ios) emin, emax,firstk,lastk 261 IF(ios/=0) THEN 262 READ(aux,*) emin, emax 263 firstk=1 264 lastk=nks 265 ENDIF 266 IF(firstk>1) kx = kx-kx(firstk) 267! 268! Since the minimum and miximum energies are given in input we can 269! sign the bands that are completely outside this range. 270! 271 is_in_range = .false. 272 DO i=1,nbnd 273 is_in_range(i) = any (e(i,1:nks) >= emin .and. e(i,1:nks) <= emax) 274 ENDDO 275! 276! Now we compute how many paths there are: nlines 277! The first point of this path: point(iline) 278! How many points are in each path: npoints(iline) 279! 280 DO n=firstk,lastk 281 IF (high_symmetry(n)) THEN 282 IF (n==1) THEN 283! 284! first point. Initialize the number of lines, and the number of point 285! and say that this line start at the first point 286! 287 nlines=1 288 npoints(1)=1 289 point(1)=1 290 ELSEIF (n==nks) THEN 291! 292! Last point. Here we save the last point of this line, but 293! do not increase the number of lines 294! 295 npoints(nlines) = npoints(nlines)+1 296 point(nlines+1)=n 297 ELSE 298! 299! Middle line. The current line has one more points, and there is a new 300! line that has to be initialized. It has one point and its first point 301! is the current k. 302! 303 npoints(nlines) = npoints(nlines)+1 304 nlines=nlines+1 305 IF (nlines>max_lines) CALL errore('plotband','too many lines',1) 306 npoints(nlines) = 1 307 point(nlines)=n 308 ENDIF 309 IF (n==1) THEN 310 WRITE( stdout,'("high-symmetry point: ",3f7.4,& 311 &" x coordinate 0.0000")') (k(i,n),i=1,3) 312 ELSE 313 WRITE( stdout,'("high-symmetry point: ",3f7.4,& 314 &" x coordinate",f9.4)') (k(i,n),i=1,3), kx(n) 315 ENDIF 316 ELSE 317! 318! This k is not an high symmetry line so we just increase the number of 319! points of this line. 320! 321 npoints(nlines) = npoints(nlines)+1 322 ENDIF 323 ENDDO 324 ! 325 WRITE(*,'("output file (gnuplot/xmgr) > ")', advance="NO") 326 READ(5,'(a)', end=25, err=25) filename 327 !save this file name for plotting projected bandst 328 filenamegnu=filename 329 IF (filename == ' ' ) THEN 330 WRITE(*,'("skipping ...")') 331 GOTO 25 332 ENDIF 333 IF (.not.exist_rap) THEN 334! 335! Here the symmetry analysis has not been done. So simply save the bands 336! on output. 337! 338 OPEN (unit=2,file=filename,form='formatted',status='unknown',& 339 iostat=ios) 340 ! draw bands 341 DO i=1,nbnd 342 IF (is_in_range(i)) THEN 343 IF (exist_proj) THEN 344 WRITE (2,'(3f10.4)') (kx(n), e(i,n), sumproj(i,n),n=firstk,lastk) 345 ELSE 346 WRITE (2,'(2f10.4)') (kx(n), e(i,n),n=firstk,lastk) 347 ENDIF 348 WRITE (2,*) 349 ENDIF 350 ENDDO 351 CLOSE (unit = 2) 352 ELSE 353! 354! In this case we write a diffent file for each line and for each 355! representation. Each file contains the bands of that representation. 356! The file is called filename.#line.#rap 357! 358! First determine for each line how many representations are there 359! in each line 360! 361 DO ilines=1,nlines 362 nrap(ilines)=0 363 DO ipoint=1,npoints(ilines)-2 364 n=point(ilines) + ipoint 365 DO ibnd=1,nbnd 366 nrap(ilines)=max(nrap(ilines),rap(ibnd,n)) 367 ENDDO 368 ENDDO 369 IF (nrap(ilines) > 12) CALL errore("plotband",& 370 "Too many representations",1) 371 ENDDO 372! 373! Then, for each line and for each representation along that line 374! 375 DO ilines=1,nlines 376 IF (nrap(ilines)==0) THEN 377! 378! Along this line the symmetry decomposition has not been done. 379! Plot all the bands as in the standard case 380! 381 filename1=trim(filename) // "." // trim(int_to_char(ilines)) 382 383 OPEN (unit=2,file=filename1,form='formatted',status='unknown',& 384 iostat=ios) 385 ! draw bands 386 DO i=1,nbnd 387 IF (is_in_range(i)) THEN 388 !!! 389 IF (exist_proj) THEN 390 WRITE (2,'(3f10.4)') (kx(n), e(i,n), sumproj(i,n), & 391 n=point(ilines),point(ilines+1)) 392 ELSE 393 !!! 394 WRITE (2,'(2f10.4)') (kx(n), e(i,n),n=point(ilines),& 395 point(ilines+1)) 396 ENDIF 397 WRITE (2,*) 398 ENDIF 399 ENDDO 400 CLOSE (unit = 2) 401 ENDIF 402 403 todo=.true. 404 DO irap=1, nrap(ilines) 405! 406! open a file 407! 408 filename1=trim(filename) // "." // trim(int_to_char(ilines)) & 409 // "." // trim(int_to_char(irap)) 410 OPEN (unit=2,file=filename1,form='formatted',status='unknown',& 411 iostat=ios) 412 IF (ios /= 0) CALL errore("plotband","opening file" & 413 //trim(filename1),1) 414! For each k point along this line selects only the bands which belong 415! to the irap representation 416 nbnd_rapk=100000 417 DO n=point(ilines)+1, point(ilines+1)-1 418 nbnd_rapk(n)=0 419 DO i=1,nbnd 420 IF (rap(i,n)==irap) THEN 421 nbnd_rapk(n) = nbnd_rapk(n) + 1 422 e_rap(nbnd_rapk(n),n)=e(i,n) 423 !!! 424 IF (exist_proj) p_rap(nbnd_rapk(n),n)=sumproj(i,n) 425 !!! 426 ENDIF 427 ENDDO 428 ENDDO 429! 430! on the two high symmetry points the representation is different. So for each 431! band choose the closest eigenvalue available. 432! 433 DO i=1,nbnd_rapk(point(ilines)+1) 434 mine=1.e8 435 DO j=1,nbnd 436 IF (abs(e_rap(i,point(ilines)+1)-e(j,point(ilines)))<mine & 437 .and. todo(j,1)) THEN 438 e_rap(i,point(ilines))=e(j,point(ilines)) 439 !!! 440 IF (exist_proj) p_rap(i,point(ilines))=sumproj(j,point(ilines)) 441 !!! 442 mine=abs( e_rap(i,point(ilines)+1)-e(j,point(ilines))) 443 jnow=j 444 ENDIF 445 ENDDO 446 todo(jnow,1)=.false. 447 ENDDO 448 DO i=1,nbnd_rapk(point(ilines+1)-1) 449 mine=1.e8 450 DO j=1,nbnd 451 IF (abs(e_rap(i,point(ilines+1)-1)- & 452 e(j,point(ilines+1)))<mine .and. todo(j,2)) THEN 453 e_rap(i,point(ilines+1))=e(j,point(ilines+1)) 454 !!! 455 IF (exist_proj) p_rap(i,point(ilines+1))=sumproj(j,point(ilines+1)) 456 !!! 457 mine=abs(e_rap(i,point(ilines+1)-1)-e(j,point(ilines+1)) ) 458 jnow=j 459 ENDIF 460 ENDDO 461 todo(jnow,2)=.false. 462 ENDDO 463 is_in_range_rap=.false. 464 DO i=1,minval(nbnd_rapk) 465 is_in_range_rap(i) = any (e_rap(i,point(ilines):point(ilines+1))& 466 >= emin .and. e(i,point(ilines):point(ilines+1)) <= emax) 467 ENDDO 468 DO i=1,minval(nbnd_rapk) 469 IF (is_in_range_rap(i)) THEN 470 !!! 471 IF (exist_proj) THEN 472 WRITE (2,'(3f10.4)') (kx(n), e_rap(i,n), p_rap(i,n), & 473 n=point(ilines),point(ilines+1)) 474 ELSE 475 !!! 476 WRITE (2,'(2f10.4)') (kx(n), e_rap(i,n), & 477 n=point(ilines),point(ilines+1)) 478 ENDIF 479 WRITE (2,*) 480 ENDIF 481 ENDDO 482 IF (minval(nbnd_rapk)==0) THEN 483 CLOSE (unit = 2,status='delete') 484 ELSE 485 CLOSE (unit = 2,status='keep') 486 ENDIF 487 ENDDO 488 ENDDO 489 490 ! if *.proj file is found, we also simply safe the data, 491 ! for the plotting of projected band. / Junfeng Qiao 492 IF (exist_proj) THEN 493 OPEN (unit=2,file=filename,form='formatted',status='unknown',& 494 iostat=ios) 495 ! draw bands 496 DO i=1,nbnd 497 IF (is_in_range(i)) THEN 498 WRITE (2,'(3f10.4)') (kx(n), e(i,n), sumproj(i,n),n=1,nks) 499 WRITE (2,*) 500 ENDIF 501 ENDDO 502 CLOSE (unit = 2) 503 ENDIF 504 ENDIF 505 WRITE(*,'("bands in gnuplot/xmgr format written to file ",a)') filename 506 ! 50725 CONTINUE 508 IF (exist_rap) THEN 509 DEALLOCATE(nbnd_rapk) 510 DEALLOCATE(e_rap) 511 DEALLOCATE(rap) 512 DEALLOCATE(k_rap) 513 DEALLOCATE(todo) 514 DEALLOCATE(is_in_range_rap) 515 ENDIF 516 IF (exist_proj) THEN 517 DEALLOCATE(sumproj) 518 DEALLOCATE(p_rap) 519 ENDIF 520 WRITE(*,'("output file (ps) > ")', advance="NO") 521 READ(5,'(a)',end=30,err=30) filename 522 IF (filename == ' ' ) THEN 523 WRITE(*,'("stopping ...")') 524 GOTO 30 525 ENDIF 526 OPEN (unit=1,file=trim(filename),form='formatted',status='unknown',& 527 iostat=ios) 528 WRITE(*,'("Efermi > ")', advance="NO") 529 READ(5,*) Ef 530 WRITE(*,'("deltaE, reference E (for tics) ")', advance="NO") 531 READ(5,*) deltaE, eref 532 ! 533 WRITE (1,'(a)') '%! PS-Adobe-1.0' 534 WRITE (1,*) '/localdict 100 dict def' 535 WRITE (1,*) 'localdict begin' 536 WRITE (1,*) '% delete next line for insertion in a LaTeX file' 537 WRITE (1,*) ' 0 0 moveto' 538 WRITE (1,*) 'gsave' 539 WRITE (1,*) '/nm {newpath moveto} def' 540 WRITE (1,*) '/riga {newpath moveto lineto stroke} def' 541 WRITE (1,*) '/banda {3 1 roll moveto {lineto} repeat stroke} def' 542 WRITE (1,*) '/dot {newpath 1 0 360 arc fill} def' 543 WRITE (1,*) '/Times-Roman findfont 12 scalefont setfont' 544 WRITE (1,*) 'currentpoint translate' 545 WRITE (1,*) '% Landscape: uncomment next line' 546 WRITE (1,*) ' 90 rotate 0 21 neg 28.451 mul translate 1.5 1.5 scale' 547 WRITE (1,*) '% Landscape: comment next line' 548 WRITE (1,*) '% 1.2 1.2 scale' 549 WRITE (1,'(2(f8.3,1x)," translate")') x0, y0 550 WRITE (1,*) '0 setgray 0.5 setlinewidth' 551 ! draw tics on axis 552 ni=nint((eref-emin)/deltaE)+1 553 nf=nint((emax-eref)/deltaE)+1 554 DO i=-ni,nf 555 etic=eref+i*deltaE 556 IF (etic >= emin .and. etic <= emax) THEN 557 WRITE (1,'(2(f8.3,1x)," moveto -5 0 rlineto stroke")') & 558 0.0,(etic-emin)*ydim/(emax-emin) 559 WRITE (1,'(2(f8.3,1x)," moveto (",f5.1,") show")') & 560 -30.,(etic-emin)*ydim/(emax-emin), etic-eref 561 ENDIF 562 ENDDO 563 ! draw the Fermi Energy 564 IF (Ef > emin .and. Ef < emax) THEN 565 WRITE (1,'("[2 4] 0 setdash newpath ",2(f8.3,1x), " moveto ")') & 566 0.0, (Ef-emin)/(emax-emin)*ydim 567 WRITE (1,'(2(f8.3,1x)," lineto stroke [] 0 setdash")') & 568 xdim, (Ef-emin)/(emax-emin)*ydim 569 ENDIF 570 ! draw axis and set clipping region 571 WRITE (1,*) '1 setlinewidth' 572 WRITE (1,'(8(f8.3,1x))') 0.0,0.0,0.0,ydim,xdim,ydim,xdim,0.0 573 WRITE (1,*) 'newpath moveto lineto lineto lineto closepath clip stroke' 574 WRITE (1,*) '0.5 setlinewidth' 575 ! draw high-symmetry lines 576 DO n=1,nks 577 IF (high_symmetry(n)) THEN 578 WRITE (1,'(4(f8.3,1x)," riga")') & 579 kx(n)*xdim/kx(nks), 0.0, kx(n)*xdim/kx(nks), ydim 580 ENDIF 581 DO i=1,nbnd 582 IF (is_in_range(i)) WRITE (1,'(2(f8.3,1x)," dot")' ) & 583 kx(n)*xdim/kx(nks), (e(i,n)-emin)*ydim/(emax-emin) 584 ENDDO 585 ENDDO 586 ! draw bands 587 ALLOCATE (k_interp(4*nks), e_interp(4*nks), coef_interp(nks,4)) 588 DO i=1,nbnd 589 IF (is_in_range(i)) THEN 590 ! No interpolation: 591 ! write (1,'(9(f8.3,1x))') ( kx(n)*xdim/kx(nks), & 592 ! (e(i,n)-emin)*ydim/(emax-emin),n=nks,1,-1) 593 ! write (1,'(i4," banda")' ) nks-1 594 ! Spline interpolation with twice as many points: 595 ! 596 ni=1 597 nf=1 598 DO nl=1,nlines 599 ni=nf 600 nf=nf + npoints(nl)-1 601 n_interp= 2*(nf-ni)+1 602 IF (n_interp < 7) CYCLE 603 DO n=1,n_interp 604 k_interp(n)=kx(ni)+(n-1)*(kx(nf)-kx(ni))/(n_interp-1) 605 ENDDO 606 DO n=ni,nf 607 e_in(n-ni+1)=e(i,n) 608 ENDDO 609 CALL spline_interpol ( kx(ni), e_in, nf-ni+1, & 610 k_interp, e_interp, n_interp ) 611 WRITE (1,'(9(f8.3,1x))') ( k_interp(n)*xdim/kx(nks), & 612 (e_interp(n)-emin)*ydim/(emax-emin),n=n_interp,1,-1) 613 WRITE (1,'(i4," banda")' ) n_interp-1 614 ENDDO 615 ENDIF 616 ENDDO 617 618 WRITE (1,*) 'grestore' 619 WRITE (1,*) '% delete next lines for insertion in a tex file' 620 WRITE (1,'(a)') '%%Page' 621 WRITE (1,*) 'showpage' 622 CLOSE (unit=1) 623 WRITE(*,'("bands in PostScript format written to file ",a)') filename 62430 CONTINUE 625 626 ! generate gnuplot script to directly draw projected band structure 627 ! Junfeng Qiao (Sep/12/2018) 628 IF (exist_proj) THEN 629 WRITE(*,'(a)', advance="NO") "output file for projected band & 630 &(gnuplot script) > " 631 READ(5,'(a)', end=35, err=35) filename 632 OPEN (unit=1,file=trim(filename),form='formatted',& 633 status='unknown',iostat=ios) 634 WRITE (1,'(a)') '#!gnuplot' 635 WRITE (1,'(a)') 'set terminal postscript portrait & 636 & enhanced color dashed lw 1 ",12"' 637 WRITE (1,'(a)') 'set output "'//trim(filename)//'_projected.ps"' 638 WRITE (1,'(a,f10.4,a)') 'set xrange [0.0:',kx(nks),']' 639 WRITE (1,'(a)') 'unset xtics' 640 WRITE (1,'(a,f12.6,a,f12.6,a)') '#set yrange [',emin,':',emax,']' 641 WRITE (1,'(a,3(f12.6,a))') 'set ytics ',emin,',',deltaE,',',emax,' ' 642 WRITE (1,'(a)') 'set ylabel "E - E_{ref} (eV)"' 643 WRITE (1,'(a)') 'set border lw 0.5' 644 WRITE (1,'(a)') "set style arrow 1 nohead front lw 0.5 lc rgb 'black'" 645 DO nl=1,nlines 646 WRITE (1,'(a,f10.4,a,f10.4,a)') 'set arrow from ',kx(point(nl)),& 647 &',graph 0 to ',kx(point(nl)),',graph 1 as 1' 648 ENDDO 649 WRITE (1,'(a)') "set title '"//trim(filename)//"_projected' noenhanced" 650 WRITE (1,'(a,f12.6,a)') & 651 &"plot '"//trim(filenamegnu)//& 652 &"' u 1:($2 - ",eref,"):3 w l palette lw 1 notitle, "//CHAR(92) 653 ! char(92) = backslash; syntax "something \" confuses the PGI compiler 654 WRITE (1,'(f12.6,a)') & 655 &Ef-eref," lt 2 lw 0.5 lc rgb 'grey50' notitle" 656 CLOSE (unit=1) 657 WRITE (*,'(a)') 'run "gnuplot '//trim(filename)//'" to get "'//& 658 &trim(filename)//'_projected.ps"' 659 WRITE (*,'(a)') 'and/or run "ps2pdf '//trim(filename)//& 660 &'_projected.ps" to get "'//trim(filename)//& 661 &'_projected.pdf"' 662 ENDIF 663 66435 CONTINUE 665 STOP 66620 WRITE(*, '("Error reading k-point # ",i4)') n 667 STOP 66822 WRITE(*, '("Error reading projection file header")') 669 STOP 67023 WRITE(*, '("Error reading projections")') 671 STOP 672 673CONTAINS 674 675SUBROUTINE spline_interpol (xin, yin, nin, xout, yout, nout) 676 677 ! xin and xout should be in increasing order, with 678 ! xout(1) <= xin(1), xout(nout) <= xin(nin) 679 680 IMPLICIT NONE 681 INTEGER, INTENT(in) :: nin, nout 682 real, INTENT(in) :: xin(nin), yin(nin), xout(nout) 683 real, INTENT(out) :: yout(nout) 684 ! work space (automatically allocated) 685 real :: d2y(nin) 686 real :: dy1, dyn 687 688 dy1 = (yin(2)-yin(1))/(xin(2)-xin(1)) 689 dyn = 0.0 690 691 CALL spline( xin, yin, nin, dy1, dyn, d2y) 692 CALL splint( nin, xin, yin, d2y, nout, xout, yout) 693 694 RETURN 695END SUBROUTINE spline_interpol 696 697SUBROUTINE spline(x, y, n, yp1, ypn, d2y) 698 699 IMPLICIT NONE 700 INTEGER, INTENT(in) :: n 701 real, INTENT(in) :: x(n), y(n), yp1, ypn 702 real, INTENT(out):: d2y(n) 703 ! work space (automatically allocated) 704 real :: work(n) 705 INTEGER :: i, k 706 real :: sig, p, qn, un 707 708 d2y(1)=-0.5 709 work(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) 710 711 DO i=2,n-1 712 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) 713 p=sig*d2y(i-1)+2.0 714 d2y(i)=(sig-1.0)/p 715 work(i)=(6.0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) & 716 /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*work(i-1))/p 717 ENDDO 718 qn=0.5 719 un=(3.0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) 720 721 d2y(n)=(un-qn*work(n-1))/(qn*d2y(n-1)+1.0) 722 DO k=n-1,1,-1 723 d2y(k)=d2y(k)*d2y(k+1)+work(k) 724 ENDDO 725 726 RETURN 727END SUBROUTINE spline 728 729 730SUBROUTINE splint (nspline, xspline, yspline, d2y, nfit, xfit, yfit) 731 732 IMPLICIT NONE 733 ! input 734 INTEGER, INTENT(in) :: nspline, nfit 735 real, INTENT(in) :: xspline(nspline), yspline(nspline), xfit(nfit), & 736 d2y(nspline) 737 real, INTENT(out) :: yfit(nfit) 738 INTEGER :: klo, khi, i 739 real :: a, b, h 740 741 IF (nspline==2) THEN 742 PRINT *, "n=",nspline,nfit 743 PRINT *, xspline 744 PRINT *, yspline 745 PRINT *, d2y 746 ENDIF 747 klo=1 748 DO i=1,nfit 749 DO khi=klo+1, nspline 750 IF(xspline(khi) >= xfit(i)) THEN 751 IF(xspline(khi-1) <= xfit(i)) THEN 752 klo = khi-1 753 ELSE 754 IF (klo == 1 .and. khi-1 == 1) THEN 755 ! the case xfit(i) < xspline(1) should not happen 756 ! but since it may be due to a numerical artifact 757 ! we just continue 758 PRINT *, ' SPLINT WARNING: xfit(i) < xspline(1)', & 759 xfit(i), xspline(1) 760 ELSE 761 STOP ' SPLINT ERROR: xfit not properly ordered' 762 ENDIF 763 ENDIF 764 h= xspline(khi) - xspline(klo) 765 a= (xspline(khi)-xfit(i))/h 766 b= (xfit(i)-xspline(klo))/h 767 768 yfit(i) = a*yspline(klo) + b*yspline(khi) & 769 + ( (a**3-a)*d2y(klo) + (b**3-b)*d2y(khi) )*h*h/6.0 770 GOTO 10 771 ENDIF 772 ENDDO 773 774 ! the case xfit(i) > xspline(nspline) should also not happen 775 ! but again it may be due to a numerical artifact 776 ! A properly chosen extrapolation formula should be used here 777 ! (and in the case xfit(i) < xspline(1) above as well) but 778 ! I am too lazy to write one - PG 779 780 PRINT *, ' SPLINT WARNING: xfit(i) > xspline(nspline)', & 781 xfit(i), xspline(nspline) 782 khi = klo+1 783 h= xspline(khi) - xspline(klo) 784 a= (xspline(khi)-xfit(i))/h 785 b= (xfit(i)-xspline(klo))/h 786 787 yfit(i) = a*yspline(klo) + b*yspline(khi) & 788 + ( (a**3-a)*d2y(klo) + (b**3-b)*d2y(khi) )*h*h/6.0 789 ! 79010 CONTINUE 791 ENDDO 792 793 RETURN 794END SUBROUTINE splint 795 796SUBROUTINE readline(aunit, inline) 797 ! read input of arbitrary length, 798 ! return a string of length at least 256. 799 ! the returning string will have at least one 800 ! whitespace, to be compatible with field_count() 801 ! Junfeng Qiao (Sep/12/2018) 802 IMPLICIT NONE 803 INTEGER, INTENT(in) :: aunit 804 CHARACTER(len=:), ALLOCATABLE, INTENT(out) :: inline 805 CHARACTER(len=:), ALLOCATABLE :: tmpline 806 INTEGER, PARAMETER :: line_buf_len=256 807 CHARACTER(len=line_buf_len) :: instr 808 LOGICAL :: set 809 INTEGER status, size 810 811 set = .true. 812 DO 813 READ(aunit,'(a)',advance='NO',iostat=status, size=size) instr 814 IF (set) THEN 815 tmpline = instr(1:size) 816 set=.false. 817 ELSE 818 tmpline = tmpline // instr(1:size) 819 ENDIF 820 IF (IS_IOSTAT_EOR(status)) exit 821 ENDDO 822 ! the inline will have at least one blank at the ending 823 IF (len_trim(tmpline) < 256) THEN 824 instr = tmpline 825 inline = instr 826 ELSE 827 inline = tmpline // ' ' 828 ENDIF 829 RETURN 830END SUBROUTINE readline 831 832END PROGRAM plotband 833 834