TABLE OF CONTENTS
ABINIT/band2eps [ Programs ]
NAME
band2eps
FUNCTION
Draws the phonon dispersion curve in Encapsuled PostScript (EPS) in black and white or in color according to the displacement participation of each atom.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (FDortu,MVeithen) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
(main routine)
OUTPUT
(main routine)
SOURCE
25 #if defined HAVE_CONFIG_H 26 #include "config.h" 27 #endif 28 29 #include "abi_common.h" 30 31 32 program band2eps 33 34 use defs_basis 35 use m_abimover 36 use m_xmpi 37 use m_abicore 38 use m_errors 39 use m_effective_potential 40 use m_multibinit_dataset 41 use m_effective_potential_file 42 use m_band2eps_dataset 43 44 use m_build_info, only : abinit_version 45 use m_io_tools, only : open_file 46 use m_fstrings, only : int2char4, tolower, inupper 47 use m_time, only : asctime 48 use m_parser, only : instrng 49 50 implicit none 51 52 !Arguments ----------------------------------- 53 !Local variables------------------------------- 54 integer,parameter :: master=0 55 character(len=fnlen) :: filnam(4) 56 real(dp) :: E,deltaE 57 integer :: comm,EmaxN,EminN,kmaxN,kminN,lastPos,lenstr,pos,posk 58 integer :: iatom,ii,imode,io,iqpt,jj,nqpt 59 integer :: nproc,my_rank 60 integer :: option,unt1,unt2,unt3 61 logical :: iam_master 62 !array 63 real(dp),allocatable :: phfrq(:),phfrqqm1(:) 64 real(dp),allocatable :: color(:,:) 65 real(dp) :: facUnit,norm,renorm 66 real(dp),allocatable :: colorAtom(:,:) 67 real(dp),allocatable :: displ(:,:) 68 type(band2eps_dataset_type) :: inp 69 character(len=500) :: message 70 character(len=strlen) :: string, raw_string 71 !scale : hold the scale for each line (dimension=nlines) 72 !qname : hold the name (gamma,R,etc..) for each extremity of line (dimension=nlines+1) 73 !nqptl : =nqpt by line (dimension=nlines) 74 !nlines : number of lines 75 !Emin is the minimum energy of the vertical axe 76 !Emax is the maximum energy of the vertical axe 77 !EminN is the minimum value of the vertical axe(in point) 78 !EmaxN is the maximum value of the vertical axe(in point) 79 !kminN is the minimum value of the horizontal axe(in point) 80 !kmaxN is the maximum value of the horizontal axe(in point) 81 !E,deltaE,pos are a work variables 82 !gradRes is the number of intervals for the graduation along vertical axe 83 84 ! ********************************************************************* 85 86 !Change communicator for I/O (mandatory!) 87 call abi_io_redirect(new_io_comm=xmpi_world) 88 89 !Initialize MPI 90 call xmpi_init() 91 comm = xmpi_world 92 93 !MPI variables 94 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 95 iam_master = (my_rank == master) 96 97 !Initialize memory profiling if it is activated 98 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0" 99 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally 100 #ifdef HAVE_MEM_PROFILING 101 call abimem_init(0) 102 #endif 103 104 !read the .file file 105 !File names refer to following files, in order: 106 !(1) Formatted input file 107 !(2) EPS graphic 108 !(3) Input phonon energies (from sortph.f) 109 !(4) Input displacements (from sortph.f) 110 write(std_out,*)' Give name for formatted input file : ' 111 read(std_in, '(a)',IOSTAT=io) filnam(1) 112 write(std_out,'(a,a)' )'- ',trim(filnam(1)) 113 114 write(std_out,*)' Give name for formatted output eps file : ' 115 read(std_in, '(a)',IOSTAT=io) filnam(2) 116 write(std_out,'(a,a)' )'- ',trim(filnam(2)) 117 118 write(std_out,*)' Give name for formatted phonon frequency file : ' 119 read(std_in, '(a)',IOSTAT=io) filnam(3) 120 write(std_out,'(a,a)' )'- ',trim(filnam(3)) 121 122 write(std_out,*)' Give name for formatted displacements file : ' 123 read(std_in, '(a)',IOSTAT=io) filnam(4) 124 write(std_out,'(a,a)' )'- ',trim(filnam(4)) 125 126 !Read the input file, and store the information in a long string of characters 127 !strlen from defs_basis module 128 write(std_out,'(a,a)') 'Opening and reading input file: ', filnam(1) 129 option=1 130 call instrng (filnam(1),lenstr,option,strlen,string, raw_string) 131 !To make case-insensitive, map characters to upper case: 132 call inupper(string(1:lenstr)) 133 134 !Read the input file 135 call invars11(inp,lenstr,string) 136 if(inp%prtout == 1) call outvars_band2eps(inp,std_out) 137 138 !Open the '.eps' file for write 139 write(std_out,'(a,a)') 'Creation of file ', filnam(2) 140 if (open_file(filnam(2),message,newunit=unt1,form="formatted",status="unknown",action="write") /= 0) then 141 ABI_ERROR(message) 142 end if 143 !Open the phonon energies file 144 if (open_file(filnam(3),message,newunit=unt2,form="formatted") /= 0) then 145 ABI_ERROR(message) 146 end if 147 if(filnam(4)/='no') then 148 ! Open the displacements file 149 if (open_file(filnam(4),message,newunit=unt3,form="formatted",status="old",action='read') /= 0) then 150 ABI_ERROR(message) 151 end if 152 end if 153 154 155 !Boundings of the plot (only the plot and not what is around) 156 EminN=6900 157 EmaxN=2400 158 kminN=2400 159 kmaxN=9600 160 161 !Allocate dynamique variables 162 ABI_MALLOC(phfrqqm1,(3*inp%natom)) 163 ABI_MALLOC(phfrq,(3*inp%natom)) 164 ABI_MALLOC(color,(3,3*inp%natom)) 165 ABI_MALLOC(colorAtom,(3,inp%natom)) 166 !colorAtom(1,1:5) : atoms contributing to red (ex : [1 0 0 0 0]) 167 !colorAtom(2,1:5) : atoms contributing to green (ex : [0 1 0 0 0]) 168 !colorAtom(3,1:5) : atoms contributing to blue (ex : [0 0 1 1 1]) 169 !tranfert color from input 170 colorAtom(1,:) = inp%red 171 colorAtom(2,:) = inp%green 172 colorAtom(3,:) = inp%blue 173 ABI_MALLOC(displ,(inp%natom,3*inp%natom)) 174 !Read end of input file 175 176 !Multiplication factor for units (from Hartree to cm-1 or THz) 177 if(inp%cunit==1) then 178 facUnit=Ha_cmm1 179 elseif(inp%cunit==2) then 180 facUnit=Ha_THz 181 else 182 end if 183 !calculate nqpt 184 nqpt=0 185 do ii=1,inp%nlines 186 nqpt=nqpt+inp%nqline(ii) 187 end do 188 !compute normalisation factor 189 renorm=0 190 do ii=1,inp%nlines 191 renorm=renorm+inp%nqline(ii)*inp%scale(ii) 192 end do 193 renorm=renorm/nqpt 194 !Calculate inp%min and inp%max 195 inp%min=inp%min/FacUnit 196 inp%max=inp%max/FacUnit 197 198 !******************************************************* 199 !Begin to write some comments in the eps file 200 !This is based to 'xfig' 201 202 write(unt1,'(a)') '% !PS-Adobe-2.0 EPSF-2.0' 203 write(unt1,'(a)') '%%Title: band.ps' 204 write(unt1,'(a)') '%%BoundingBox: 0 0 581 310' 205 write(unt1,'(a)') '%%Magnification: 1.0000' 206 207 write(unt1,'(a)') '/$F2psDict 200 dict def' 208 write(unt1,'(a)') '$F2psDict begin' 209 write(unt1,'(a)') '$F2psDict /mtrx matrix put' 210 write(unt1,'(a)') '/col-1 {0 setgray} bind def' 211 write(unt1,'(a)') '/col0 {0.000 0.000 0.000 srgb} bind def' 212 write(unt1,'(a)') 'end' 213 write(unt1,'(a)') 'save' 214 write(unt1,'(a)') 'newpath 0 310 moveto 0 0 lineto 581 0 lineto 581 310 lineto closepath clip newpath' 215 write(unt1,'(a)') '-36.0 446.0 translate' 216 write(unt1,'(a)') '1 -1 scale' 217 218 write(unt1,'(a)') '/cp {closepath} bind def' 219 write(unt1,'(a)') '/ef {eofill} bind def' 220 write(unt1,'(a)') '/gr {grestore} bind def' 221 write(unt1,'(a)') '/gs {gsave} bind def' 222 write(unt1,'(a)') '/sa {save} bind def' 223 write(unt1,'(a)') '/rs {restore} bind def' 224 write(unt1,'(a)') '/l {lineto} bind def' 225 write(unt1,'(a)') '/m {moveto} bind def' 226 write(unt1,'(a)') '/rm {rmoveto} bind def' 227 write(unt1,'(a)') '/n {newpath} bind def' 228 write(unt1,'(a)') '/s {stroke} bind def' 229 write(unt1,'(a)') '/sh {show} bind def' 230 write(unt1,'(a)') '/slc {setlinecap} bind def' 231 write(unt1,'(a)') '/slj {setlinejoin} bind def' 232 write(unt1,'(a)') '/slw {setlinewidth} bind def' 233 write(unt1,'(a)') '/srgb {setrgbcolor} bind def' 234 write(unt1,'(a)') '/rot {rotate} bind def' 235 write(unt1,'(a)') '/sc {scale} bind def' 236 write(unt1,'(a)') '/sd {setdash} bind def' 237 write(unt1,'(a)') '/ff {findfont} bind def' 238 write(unt1,'(a)') '/sf {setfont} bind def' 239 write(unt1,'(a)') '/scf {scalefont} bind def' 240 write(unt1,'(a)') '/sw {stringwidth} bind def' 241 write(unt1,'(a)') '/tr {translate} bind def' 242 write(unt1,'(a)') '/tnt {dup dup currentrgbcolor' 243 244 write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add' 245 write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add' 246 write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add srgb}' 247 write(unt1,'(a)') 'bind def' 248 write(unt1,'(a)') '/shd {dup dup currentrgbcolor 4 -2 roll mul 4 -2 roll mul' 249 write(unt1,'(a)') ' 4 -2 roll mul srgb} bind def' 250 write(unt1,'(a)') '/$F2psBegin {$F2psDict begin /$F2psEnteredState save def} def' 251 write(unt1,'(a)') '/$F2psEnd {$F2psEnteredState restore end} def' 252 write(unt1,'(a)') '$F2psBegin' 253 write(unt1,'(a)') '%%Page: 1 1' 254 write(unt1,'(a)') '10 setmiterlimit' 255 write(unt1,'(a)') '0.06000 0.06000 sc' 256 257 !**************************************************************** 258 !Begin of the intelligible part of the postcript document 259 260 write(unt1,'(a)') '%**************************************' 261 !**************************************************************** 262 !Draw the box containing the plot 263 write(unt1,'(a)') '%****Big Box****' 264 write(unt1,'(a)') '12 slw' 265 write(unt1,'(a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ', EmaxN,& 266 & ' m ', kmaxN,' ', EmaxN, ' l ', & 267 & kmaxN,' ', EminN, ' l ', kminN,' ', EminN, ' l' 268 write(unt1,'(a)') 'cp gs col0 s gr' 269 270 !**************************************************************** 271 !Write unit on the middle left of the vertical axe 272 write(unt1,'(a)') '%****Units****' 273 274 if(inp%cunit==1) then 275 ! 1/lambda 276 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 277 write(unt1,'(a)') '1425 5650 m' 278 write(unt1,'(3a)') 'gs 1 -1 sc 90.0 rot (Frequency ',achar(92),'(cm) col0 sh gr' 279 ! cm-1 280 write(unt1,'(a)') '/Times-Roman ff 200.00 scf sf' 281 write(unt1,'(a)') '1325 4030 m' 282 write(unt1,'(a)') 'gs 1 -1 sc 90.0 rot (-1) col0 sh gr' 283 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 284 write(unt1,'(a)') '1425 3850 m' 285 write(unt1,'(3a)') 'gs 1 -1 sc 90.0 rot (',achar(92),')) col0 sh gr' 286 else 287 ! Freq 288 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 289 write(unt1,'(a)') '825 4850 m' 290 write(unt1,'(a)') 'gs 1 -1 sc 90.0 rot (Freq) col0 sh gr' 291 ! THz 292 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 293 write(unt1,'(a)') '825 4350 m' 294 write(unt1,'(a)') 'gs 1 -1 sc 90.0 rot (THz) col0 sh gr' 295 end if 296 !***************************************************************** 297 !Write graduation on the vertical axe 298 write(unt1,'(a)') '%****Vertical graduation****' 299 deltaE=(inp%max-inp%min)/inp%ngrad 300 301 E=inp%min 302 do 303 ! do E=inp%min,(inp%max-deltaE/2),deltaE 304 if (E >= (inp%max-deltaE/2)-tol6) exit 305 pos=int(((EminN-EmaxN)*E & 306 & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max)) 307 308 ! write the value of energy(or frequence) 309 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 310 write(unt1,'(i4,a,i4,a)') kminN-800,' ',pos+60,' m' !-1300 must be CHANGED 311 ! as a function of the width of E 312 write(unt1,'(a,i6,a)') 'gs 1 -1 sc (', nint(E*facUnit),') col0 sh gr' 313 314 ! write a little bar 315 write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kminN+100,' ', pos, ' l' 316 write(unt1,'(a)') 'gs col0 s gr ' 317 318 E = E+deltaE 319 end do 320 321 !do the same thing for E=inp%max (floating point error) 322 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 323 write(unt1,'(i4,a,i4,a)') kminN-800,' ',EmaxN+60,' m' !-1300 must be changed as E 324 write(unt1,'(a,i6,a)') 'gs 1 -1 sc (', nint(inp%max*facUnit),') col0 sh gr' 325 326 327 !draw zero line 328 E=0 329 pos=int(((EminN-EmaxN)*E & 330 & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max)) 331 write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kmaxN,' ', pos, ' l' 332 write(unt1,'(a)') 'gs col0 s gr ' 333 334 335 !****************************************************** 336 !draw legend of horizontal axe 337 !+vertical line 338 339 write(unt1,'(a)') '%****Horizontal graduation****' 340 341 lastPos=kminN 342 343 do ii=0,inp%nlines 344 345 if(ii/=0) then 346 posk=int(((kminN-kmaxN)*(inp%nqline(ii))) & 347 & *inp%scale(ii)/renorm/(-nqpt)) 348 else 349 posk=0 350 end if 351 352 posk=posk+lastPos 353 lastPos=posk 354 355 if(tolower(inp%qpoint_name(ii+1))=='gamma') then !GAMMA 356 write(unt1,'(a)') '/Symbol ff 270.00 scf sf' 357 write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m' 358 write(unt1,'(a)') 'gs 1 -1 sc (G) col0 sh gr' 359 elseif(tolower(inp%qpoint_name(ii+1))=='lambda') then !LAMBDA 360 write(unt1,'(a)') '/Symbol ff 270.00 scf sf' 361 write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m' 362 write(unt1,'(a)') 'gs 1 -1 sc (L) col0 sh gr' 363 else !autre 364 write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf' 365 write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m' 366 write(unt1,'(a,a1,a)') 'gs 1 -1 sc (',inp%qpoint_name(ii+1),') col0 sh gr' 367 end if 368 369 ! draw vertical line 370 write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', posk,' ',EminN ,' m ', posk,' ', EmaxN, ' l' 371 write(unt1,'(a)') 'gs col0 s gr ' 372 373 end do 374 375 !*********************************************************** 376 !Write the bands (the most important part actually) 377 378 write(unt1,'(a)') '%****Write Bands****' 379 380 lastPos=kminN 381 382 read(unt2,*) (phfrqqm1(ii),ii=1,3*inp%natom) 383 384 do jj=1,inp%nlines 385 do iqpt=1,inp%nqline(jj) 386 read(unt2,*) (phfrq(ii),ii=1,3*inp%natom) 387 do imode=1,3*inp%natom 388 389 if(filnam(4)/='no') then !calculate the color else in black and white 390 do iatom=1,inp%natom 391 read(unt3,*) displ(iatom,imode) 392 end do 393 ! normalize displ 394 norm=0 395 do iatom=1,inp%natom 396 norm=norm+displ(iatom,imode) 397 end do 398 399 do iatom=1,inp%natom 400 displ(iatom,imode)=displ(iatom,imode)/norm 401 end do 402 403 ! Treat color 404 color(:,imode)=0 405 do ii=1,inp%natom 406 ! Red 407 color(1,imode)=color(1,imode)+displ(ii,imode)*colorAtom(1,ii) 408 ! Green 409 color(2,imode)=color(2,imode)+displ(ii,imode)*colorAtom(2,ii) 410 ! Blue 411 color(3,imode)=color(3,imode)+displ(ii,imode)*colorAtom(3,ii) 412 end do 413 end if 414 415 pos=int(((EminN-EmaxN)*phfrqqm1(imode) & 416 & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max)) 417 418 posk=int(((kminN-kmaxN)*(iqpt-1) & 419 & *inp%scale(jj)/renorm/(-nqpt))) 420 posk=posk+lastPos 421 422 write(unt1,'(a,i5,a,i5,a)') 'n ',posk,' ',pos,' m' 423 424 pos=int(((EminN-EmaxN)*phfrq(imode) & 425 & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max)) 426 posk=int(((kminN-kmaxN)*(iqpt) & 427 & *inp%scale(jj)/renorm/(-nqpt))) 428 posk=posk+lastPos 429 write(unt1,'(i5,a,i5,a)') posk,' ',pos,' l gs' 430 431 432 if(filnam(4)/='no') then !(in color) 433 write(unt1,'(f6.3,a,f6.3,a,f6.3,a)') color(1,imode),' ', & 434 & color(2,imode),' ',color(3,imode), ' srgb s gr' 435 else 436 write(unt1,'(f6.3,a,f6.3,a,f6.3,a)') 0.0,' ', & 437 & 0.0,' ',0.0, ' srgb s gr' 438 end if 439 440 441 end do 442 443 phfrqqm1=phfrq 444 445 end do 446 lastPos=posk 447 448 end do 449 450 451 !********************************************************** 452 !Ending the poscript document 453 write(unt1,'(a)') '$F2psEnd' 454 write(unt1,'(a)') 'rs' 455 456 !call abinit_doctor("__band2eps") 457 458 call band2eps_dtset_free(inp) 459 460 call xmpi_end() 461 462 end program band2eps