Line data Source code
1 : !!****p* ABINIT/band2eps
2 : !! NAME
3 : !! band2eps
4 : !!
5 : !! FUNCTION
6 : !! Draws the phonon dispersion curve in Encapsuled PostScript (EPS)
7 : !! in black and white or in color according to the displacement participation
8 : !! of each atom.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 1999-2026 ABINIT group (FDortu,MVeithen)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! INPUTS
18 : !! (main routine)
19 : !!
20 : !! OUTPUT
21 : !! (main routine)
22 : !!
23 : !! SOURCE
24 :
25 : #if defined HAVE_CONFIG_H
26 : #include "config.h"
27 : #endif
28 :
29 : #include "abi_common.h"
30 :
31 :
32 4 : program band2eps
33 :
34 4 : 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 4 : real(dp),allocatable :: phfrq(:),phfrqqm1(:)
64 4 : real(dp),allocatable :: color(:,:)
65 : real(dp) :: facUnit,norm,renorm
66 4 : real(dp),allocatable :: colorAtom(:,:)
67 4 : real(dp),allocatable :: displ(:,:)
68 4 : 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 4 : call abi_io_redirect(new_io_comm=xmpi_world)
88 :
89 : !Initialize MPI
90 4 : call xmpi_init()
91 4 : comm = xmpi_world
92 :
93 : !MPI variables
94 4 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
95 4 : 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 4 : write(std_out,*)' Give name for formatted input file : '
111 4 : read(std_in, '(a)',IOSTAT=io) filnam(1)
112 4 : write(std_out,'(a,a)' )'- ',trim(filnam(1))
113 :
114 4 : write(std_out,*)' Give name for formatted output eps file : '
115 4 : read(std_in, '(a)',IOSTAT=io) filnam(2)
116 4 : write(std_out,'(a,a)' )'- ',trim(filnam(2))
117 :
118 4 : write(std_out,*)' Give name for formatted phonon frequency file : '
119 4 : read(std_in, '(a)',IOSTAT=io) filnam(3)
120 4 : write(std_out,'(a,a)' )'- ',trim(filnam(3))
121 :
122 4 : write(std_out,*)' Give name for formatted displacements file : '
123 4 : read(std_in, '(a)',IOSTAT=io) filnam(4)
124 4 : 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 4 : write(std_out,'(a,a)') 'Opening and reading input file: ', filnam(1)
129 4 : option=1
130 4 : call instrng (filnam(1),lenstr,option,strlen,string, raw_string)
131 : !To make case-insensitive, map characters to upper case:
132 4 : call inupper(string(1:lenstr))
133 :
134 : !Read the input file
135 4 : call invars11(inp,lenstr,string)
136 4 : if(inp%prtout == 1) call outvars_band2eps(inp,std_out)
137 :
138 : !Open the '.eps' file for write
139 4 : write(std_out,'(a,a)') 'Creation of file ', filnam(2)
140 4 : if (open_file(filnam(2),message,newunit=unt1,form="formatted",status="unknown",action="write") /= 0) then
141 0 : ABI_ERROR(message)
142 : end if
143 : !Open the phonon energies file
144 4 : if (open_file(filnam(3),message,newunit=unt2,form="formatted") /= 0) then
145 0 : ABI_ERROR(message)
146 : end if
147 4 : if(filnam(4)/='no') then
148 : ! Open the displacements file
149 2 : if (open_file(filnam(4),message,newunit=unt3,form="formatted",status="old",action='read') /= 0) then
150 0 : 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 4 : EminN=6900
157 4 : EmaxN=2400
158 4 : kminN=2400
159 4 : kmaxN=9600
160 :
161 : !Allocate dynamique variables
162 4 : ABI_MALLOC(phfrqqm1,(3*inp%natom))
163 4 : ABI_MALLOC(phfrq,(3*inp%natom))
164 4 : ABI_MALLOC(color,(3,3*inp%natom))
165 4 : 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 18 : colorAtom(1,:) = inp%red
171 18 : colorAtom(2,:) = inp%green
172 18 : colorAtom(3,:) = inp%blue
173 4 : 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 4 : if(inp%cunit==1) then
178 4 : facUnit=Ha_cmm1
179 0 : elseif(inp%cunit==2) then
180 0 : facUnit=Ha_THz
181 : else
182 : end if
183 : !calculate nqpt
184 4 : nqpt=0
185 26 : do ii=1,inp%nlines
186 26 : nqpt=nqpt+inp%nqline(ii)
187 : end do
188 : !compute normalisation factor
189 4 : renorm=0
190 26 : do ii=1,inp%nlines
191 26 : renorm=renorm+inp%nqline(ii)*inp%scale(ii)
192 : end do
193 4 : renorm=renorm/nqpt
194 : !Calculate inp%min and inp%max
195 4 : inp%min=inp%min/FacUnit
196 4 : 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 4 : write(unt1,'(a)') '% !PS-Adobe-2.0 EPSF-2.0'
203 4 : write(unt1,'(a)') '%%Title: band.ps'
204 4 : write(unt1,'(a)') '%%BoundingBox: 0 0 581 310'
205 4 : write(unt1,'(a)') '%%Magnification: 1.0000'
206 :
207 4 : write(unt1,'(a)') '/$F2psDict 200 dict def'
208 4 : write(unt1,'(a)') '$F2psDict begin'
209 4 : write(unt1,'(a)') '$F2psDict /mtrx matrix put'
210 4 : write(unt1,'(a)') '/col-1 {0 setgray} bind def'
211 4 : write(unt1,'(a)') '/col0 {0.000 0.000 0.000 srgb} bind def'
212 4 : write(unt1,'(a)') 'end'
213 4 : write(unt1,'(a)') 'save'
214 4 : write(unt1,'(a)') 'newpath 0 310 moveto 0 0 lineto 581 0 lineto 581 310 lineto closepath clip newpath'
215 4 : write(unt1,'(a)') '-36.0 446.0 translate'
216 4 : write(unt1,'(a)') '1 -1 scale'
217 :
218 4 : write(unt1,'(a)') '/cp {closepath} bind def'
219 4 : write(unt1,'(a)') '/ef {eofill} bind def'
220 4 : write(unt1,'(a)') '/gr {grestore} bind def'
221 4 : write(unt1,'(a)') '/gs {gsave} bind def'
222 4 : write(unt1,'(a)') '/sa {save} bind def'
223 4 : write(unt1,'(a)') '/rs {restore} bind def'
224 4 : write(unt1,'(a)') '/l {lineto} bind def'
225 4 : write(unt1,'(a)') '/m {moveto} bind def'
226 4 : write(unt1,'(a)') '/rm {rmoveto} bind def'
227 4 : write(unt1,'(a)') '/n {newpath} bind def'
228 4 : write(unt1,'(a)') '/s {stroke} bind def'
229 4 : write(unt1,'(a)') '/sh {show} bind def'
230 4 : write(unt1,'(a)') '/slc {setlinecap} bind def'
231 4 : write(unt1,'(a)') '/slj {setlinejoin} bind def'
232 4 : write(unt1,'(a)') '/slw {setlinewidth} bind def'
233 4 : write(unt1,'(a)') '/srgb {setrgbcolor} bind def'
234 4 : write(unt1,'(a)') '/rot {rotate} bind def'
235 4 : write(unt1,'(a)') '/sc {scale} bind def'
236 4 : write(unt1,'(a)') '/sd {setdash} bind def'
237 4 : write(unt1,'(a)') '/ff {findfont} bind def'
238 4 : write(unt1,'(a)') '/sf {setfont} bind def'
239 4 : write(unt1,'(a)') '/scf {scalefont} bind def'
240 4 : write(unt1,'(a)') '/sw {stringwidth} bind def'
241 4 : write(unt1,'(a)') '/tr {translate} bind def'
242 4 : write(unt1,'(a)') '/tnt {dup dup currentrgbcolor'
243 :
244 4 : write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add'
245 4 : write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add'
246 4 : write(unt1,'(a)') '4 -2 roll dup 1 exch sub 3 -1 roll mul add srgb}'
247 4 : write(unt1,'(a)') 'bind def'
248 4 : write(unt1,'(a)') '/shd {dup dup currentrgbcolor 4 -2 roll mul 4 -2 roll mul'
249 4 : write(unt1,'(a)') ' 4 -2 roll mul srgb} bind def'
250 4 : write(unt1,'(a)') '/$F2psBegin {$F2psDict begin /$F2psEnteredState save def} def'
251 4 : write(unt1,'(a)') '/$F2psEnd {$F2psEnteredState restore end} def'
252 4 : write(unt1,'(a)') '$F2psBegin'
253 4 : write(unt1,'(a)') '%%Page: 1 1'
254 4 : write(unt1,'(a)') '10 setmiterlimit'
255 4 : write(unt1,'(a)') '0.06000 0.06000 sc'
256 :
257 : !****************************************************************
258 : !Begin of the intelligible part of the postcript document
259 :
260 4 : write(unt1,'(a)') '%**************************************'
261 : !****************************************************************
262 : !Draw the box containing the plot
263 4 : write(unt1,'(a)') '%****Big Box****'
264 4 : write(unt1,'(a)') '12 slw'
265 4 : write(unt1,'(a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ', EmaxN,&
266 4 : & ' m ', kmaxN,' ', EmaxN, ' l ', &
267 8 : & kmaxN,' ', EminN, ' l ', kminN,' ', EminN, ' l'
268 4 : write(unt1,'(a)') 'cp gs col0 s gr'
269 :
270 : !****************************************************************
271 : !Write unit on the middle left of the vertical axe
272 4 : write(unt1,'(a)') '%****Units****'
273 :
274 4 : if(inp%cunit==1) then
275 : ! 1/lambda
276 4 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
277 4 : write(unt1,'(a)') '1425 5650 m'
278 4 : write(unt1,'(3a)') 'gs 1 -1 sc 90.0 rot (Frequency ',achar(92),'(cm) col0 sh gr'
279 : ! cm-1
280 4 : write(unt1,'(a)') '/Times-Roman ff 200.00 scf sf'
281 4 : write(unt1,'(a)') '1325 4030 m'
282 4 : write(unt1,'(a)') 'gs 1 -1 sc 90.0 rot (-1) col0 sh gr'
283 4 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
284 4 : write(unt1,'(a)') '1425 3850 m'
285 4 : write(unt1,'(3a)') 'gs 1 -1 sc 90.0 rot (',achar(92),')) col0 sh gr'
286 : else
287 : ! Freq
288 0 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
289 0 : write(unt1,'(a)') '825 4850 m'
290 0 : write(unt1,'(a)') 'gs 1 -1 sc 90.0 rot (Freq) col0 sh gr'
291 : ! THz
292 0 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
293 0 : write(unt1,'(a)') '825 4350 m'
294 0 : 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 4 : write(unt1,'(a)') '%****Vertical graduation****'
299 4 : deltaE=(inp%max-inp%min)/inp%ngrad
300 :
301 4 : E=inp%min
302 24 : do
303 : ! do E=inp%min,(inp%max-deltaE/2),deltaE
304 28 : if (E >= (inp%max-deltaE/2)-tol6) exit
305 : pos=int(((EminN-EmaxN)*E &
306 24 : & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max))
307 :
308 : ! write the value of energy(or frequence)
309 24 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
310 24 : 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 24 : write(unt1,'(a,i6,a)') 'gs 1 -1 sc (', nint(E*facUnit),') col0 sh gr'
313 :
314 : ! write a little bar
315 24 : write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kminN+100,' ', pos, ' l'
316 24 : write(unt1,'(a)') 'gs col0 s gr '
317 :
318 24 : E = E+deltaE
319 : end do
320 :
321 : !do the same thing for E=inp%max (floating point error)
322 4 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
323 4 : write(unt1,'(i4,a,i4,a)') kminN-800,' ',EmaxN+60,' m' !-1300 must be changed as E
324 4 : write(unt1,'(a,i6,a)') 'gs 1 -1 sc (', nint(inp%max*facUnit),') col0 sh gr'
325 :
326 :
327 : !draw zero line
328 4 : E=0
329 : pos=int(((EminN-EmaxN)*E &
330 4 : & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max))
331 4 : write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', kminN,' ',pos ,' m ', kmaxN,' ', pos, ' l'
332 4 : write(unt1,'(a)') 'gs col0 s gr '
333 :
334 :
335 : !******************************************************
336 : !draw legend of horizontal axe
337 : !+vertical line
338 :
339 4 : write(unt1,'(a)') '%****Horizontal graduation****'
340 :
341 4 : lastPos=kminN
342 :
343 30 : do ii=0,inp%nlines
344 :
345 26 : if(ii/=0) then
346 : posk=int(((kminN-kmaxN)*(inp%nqline(ii))) &
347 22 : & *inp%scale(ii)/renorm/(-nqpt))
348 : else
349 4 : posk=0
350 : end if
351 :
352 26 : posk=posk+lastPos
353 26 : lastPos=posk
354 :
355 26 : if(tolower(inp%qpoint_name(ii+1))=='gamma') then !GAMMA
356 8 : write(unt1,'(a)') '/Symbol ff 270.00 scf sf'
357 8 : write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
358 8 : write(unt1,'(a)') 'gs 1 -1 sc (G) col0 sh gr'
359 18 : elseif(tolower(inp%qpoint_name(ii+1))=='lambda') then !LAMBDA
360 0 : write(unt1,'(a)') '/Symbol ff 270.00 scf sf'
361 0 : write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
362 0 : write(unt1,'(a)') 'gs 1 -1 sc (L) col0 sh gr'
363 : else !autre
364 18 : write(unt1,'(a)') '/Times-Roman ff 270.00 scf sf'
365 18 : write(unt1,'(i4,a,i4,a)') posk-100,' ', 7150, ' m'
366 62 : 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 26 : write(unt1,'(a,i4,a,i4,a,i4,a,i4,a)') 'n ', posk,' ',EminN ,' m ', posk,' ', EmaxN, ' l'
371 30 : write(unt1,'(a)') 'gs col0 s gr '
372 :
373 : end do
374 :
375 : !***********************************************************
376 : !Write the bands (the most important part actually)
377 :
378 4 : write(unt1,'(a)') '%****Write Bands****'
379 :
380 4 : lastPos=kminN
381 :
382 46 : read(unt2,*) (phfrqqm1(ii),ii=1,3*inp%natom)
383 :
384 26 : do jj=1,inp%nlines
385 242 : do iqpt=1,inp%nqline(jj)
386 2260 : read(unt2,*) (phfrq(ii),ii=1,3*inp%natom)
387 2260 : do imode=1,3*inp%natom
388 :
389 2040 : if(filnam(4)/='no') then !calculate the color else in black and white
390 7200 : do iatom=1,inp%natom
391 7200 : read(unt3,*) displ(iatom,imode)
392 : end do
393 : ! normalize displ
394 1200 : norm=0
395 7200 : do iatom=1,inp%natom
396 7200 : norm=norm+displ(iatom,imode)
397 : end do
398 :
399 7200 : do iatom=1,inp%natom
400 7200 : displ(iatom,imode)=displ(iatom,imode)/norm
401 : end do
402 :
403 : ! Treat color
404 4800 : color(:,imode)=0
405 7200 : do ii=1,inp%natom
406 : ! Red
407 6000 : color(1,imode)=color(1,imode)+displ(ii,imode)*colorAtom(1,ii)
408 : ! Green
409 6000 : color(2,imode)=color(2,imode)+displ(ii,imode)*colorAtom(2,ii)
410 : ! Blue
411 7200 : 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 2040 : & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max))
417 :
418 : posk=int(((kminN-kmaxN)*(iqpt-1) &
419 2040 : & *inp%scale(jj)/renorm/(-nqpt)))
420 2040 : posk=posk+lastPos
421 :
422 2040 : write(unt1,'(a,i5,a,i5,a)') 'n ',posk,' ',pos,' m'
423 :
424 : pos=int(((EminN-EmaxN)*phfrq(imode) &
425 2040 : & +EmaxN*inp%min -EminN*inp%max)/(inp%min-inp%max))
426 : posk=int(((kminN-kmaxN)*(iqpt) &
427 2040 : & *inp%scale(jj)/renorm/(-nqpt)))
428 2040 : posk=posk+lastPos
429 2040 : write(unt1,'(i5,a,i5,a)') posk,' ',pos,' l gs'
430 :
431 :
432 2260 : if(filnam(4)/='no') then !(in color)
433 1200 : write(unt1,'(f6.3,a,f6.3,a,f6.3,a)') color(1,imode),' ', &
434 2400 : & color(2,imode),' ',color(3,imode), ' srgb s gr'
435 : else
436 840 : write(unt1,'(f6.3,a,f6.3,a,f6.3,a)') 0.0,' ', &
437 1680 : & 0.0,' ',0.0, ' srgb s gr'
438 : end if
439 :
440 :
441 : end do
442 :
443 2502 : phfrqqm1=phfrq
444 :
445 : end do
446 26 : lastPos=posk
447 :
448 : end do
449 :
450 :
451 : !**********************************************************
452 : !Ending the poscript document
453 4 : write(unt1,'(a)') '$F2psEnd'
454 4 : write(unt1,'(a)') 'rs'
455 :
456 : !call abinit_doctor("__band2eps")
457 :
458 4 : call band2eps_dtset_free(inp)
459 :
460 4 : call xmpi_end()
461 :
462 0 : end program band2eps
463 : !!***
|
