Dear QE Users, I am using Espresso-5.2.0 in Ubuntu 14.04 and facing an issue about plotting the post-processing calculation.
I have obtained the charge density files to be plotted. However, when I tried to run the command of plotrho.x, this following error appeared: At line 44 of file plotrho.f90 (unit = 1, file = 'MAPIss.rho.dat') Fortran runtime error: Bad real number in item 13 of list input I had checked the plotrho.f90 and could not identify a way to get rid of the problem. FYI, I have the output files of pp.x calculation without having a problem neither with memory usage or crash during the running process. Please kindly find both the input and plotrho.f90 files. Input file: &INPUTPP outdir='/home/efidwiindari/MAPI/', prefix='MAPIss' plot_num=0, filplot='MAPIss.rho.dat', / &PLOT nfile=1, iflag=3, output_format=5, fileout='MAPIsscontour3D.xsf' e1(1)=1.0, e1(2)=0.0, e1(3)=0.0, e2(1)=0.0, e2(2)=1.0, e2(3)=0.0, e3(1)=0.0, e3(2)=0.0, e3(3)=1.0, x0(1)=0.0, x0(2)=0.0, x0(3)=0.0, nx=45, ny=45, nz=100, / plotrho.f90 file is attached below as the file is quite long. PS : line number 44 is about unit vector which I do not have any idea why the problem arose. Furthermore, I did three calculations for three different systems and unfortunately the problems were still the same. Any comment and/or suggestions would be really appreciated. Many thanks in advance. Best Regards, Efi Dwi Indari Research Assistant, Bandung Institute of Technology
! ! Copyright (C) 2001-2007 Quantum ESPRESSO group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! ! !----------------------------------------------------------------------- PROGRAM plotrho !----------------------------------------------------------------------- ! 2D contour plot - logarithmically or linearly spaced levels ! - Postscript printable output ! if " cplot" is called: ! - contour lines plus gray levels ! - negative values are shaded ! if "psplot" is called: ! - contour lines of various kinds (solid, dashed, etc) ! IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER, PARAMETER :: stdout=6 ! for spline interpolation using essl toutines INTEGER, PARAMETER :: nwrk = 10000 real(DP) :: wrk (nwrk) INTEGER, ALLOCATABLE :: ityp (:) INTEGER :: nxi, nyi, nx, ny, i, j, k, nlevels, na, nat real(DP), ALLOCATABLE :: rhoi(:,:), xi(:), yi(:) real(DP), ALLOCATABLE :: rhoo(:,:), x (:), y (:) real(DP), ALLOCATABLE :: z (:) real(DP) :: xmin, xmax, ymin, ymax, rhomin, rhomax, rhoomin, rhoomax real(DP) :: xdim, ydim, xs, ys, r0 (3), tau1 (3), tau2 (3) real(DP), ALLOCATABLE :: tau (:,:) real(DP) :: at (3, 3), a0 CHARACTER (len=256) :: filename, fileout, ans * 1 LOGICAL :: logarithmic_scale CALL get_file ( filename ) OPEN (unit = 1, file = filename, form = 'formatted', status = 'old') READ (1, * ) nxi, nyi ALLOCATE ( xi(0:nxi), yi(0:nyi), rhoi(0:nxi,0:nyi) ) READ (1, * ) (xi (i), i = 0, nxi) READ (1, * ) (yi (j), j = 0, nyi) READ (1, * ) ( (rhoi (i, j), i = 0, nxi), j = 0, nyi) READ (1, * ) r0 READ (1, * ) tau1 READ (1, * ) tau2 READ (1, * ) nat IF (nat < 0 .or. nat > 1000000) THEN WRITE( stdout, '("Error: unlikely number of atoms ",i4)') nat STOP ENDIF ALLOCATE (tau (3,nat), ityp(nat) ) READ (1, * ) ( (tau (j, na), j = 1, 3), ityp (na), na = 1, nat) READ (1, * ) a0 READ (1, * ) at CLOSE (unit = 1) ! WRITE( stdout, '("r0 : ",3f8.4)') r0 WRITE( stdout, '("tau1 : ",3f8.4)') tau1 WRITE( stdout, '("tau2 : ",3f8.4)') tau2 ! WRITE( stdout, '("read",i4," atomic positions")') nat ! WRITE( stdout,'("Atomic positions:")') ! WRITE( stdout,'(3f8.4)') ( (tau(j,na),j=1,3),na=1,nat) WRITE( stdout, '("output file > ",$)') READ (5, '(a)') fileout WRITE( stdout, '("Read ",i3," *",i3," grid")') nxi+1, nyi+1 #ifdef __ESSL ! ! interpolation implemented only for ESSL routines ... ! WRITE( stdout, '("nx, ny (output) > ",$)') READ (5, * ) nx, ny #else nx = nxi ny = nyi #endif ALLOCATE ( x(0:nx), y(0:ny), rhoo(0:nx,0:ny) ) xmin = xi (0) xmax = xi (nxi) DO i = 0, nx x (i) = (xi (nxi) - xi (0) ) * dble (i) / dble (nx) ENDDO ymin = yi (0) ymax = yi (nyi) DO i = 0, ny y (i) = (yi (nyi) - yi (0) ) * dble (i) / dble (ny) ENDDO #ifdef __ESSL CALL dcsin2 (xi, yi, rhoi, nxi + 1, nyi + 1, nxi + 1, x, y, & nx + 1, ny + 1, rhoo, nx + 1, wrk, nwrk) #else rhoo (0:nx, 0:ny) = rhoi (0:nx, 0:ny) #endif rhomin = minval (rhoo(0:nx, 0:ny)) rhomax = maxval (rhoo(0:nx, 0:ny)) IF (rhomin > 0.d0) THEN WRITE( stdout,'("Logarithmic scale (y/n)? > ",$)') READ (5, '(a)') ans logarithmic_scale = ans/='n'.and.ans/='N' ELSE logarithmic_scale = .false. ENDIF 10 CONTINUE WRITE( stdout, '("Bounds: ",2f12.6)') rhomin, rhomax WRITE( stdout, '("min, max, # of levels > ",$)') READ (5, * ) rhoomin, rhoomax, nlevels IF ( rhoomax <= rhoomin .or. & rhoomin >= rhomax .or. rhoomax <= rhomin ) THEN WRITE( stdout, '("Out of Bounds! try again")') GOTO 10 ENDIF IF (nlevels > 1000) THEN WRITE( stdout, '("Are you sure you really need ",i8," levels?")') & nlevels ELSEIF (nlevels < 1) THEN WRITE( stdout, '("Too few levels! assuming 1 level")') nlevels = 1 ENDIF ALLOCATE (z(0:nlevels)) IF (logarithmic_scale) THEN DO k = 0, nlevels - 1 z (k) = exp (log (rhoomin) + (log (rhoomax) - log (rhoomin) ) & * dble (k) / (nlevels - 1) ) ENDDO ELSE DO k = 0, nlevels - 1 z (k) = rhoomin + (rhoomax - rhoomin) * dble (k) / (nlevels - 1) ENDDO ENDIF z (nlevels) = z (nlevels - 1) xdim = 15.0d0 * (xmax - xmin) / sqrt ( (xmax - xmin) **2 + (ymax - ymin) **2) ydim = 15.0d0 * (ymax - ymin) / sqrt ( (xmax - xmin) **2 + (ymax - ymin) **2) xs = 4.0d0 ys = 3.0d0 ! uncomment the call to "cplot" if you want contour lines, ! plus gray levels and shading for negative values CALL cplot (rhoo, nx, ny, x, xmin, xmax, nx, y, ymin, ymax, & ny, nlevels, z, xdim, ydim, xs, ys, filename, fileout) ! uncomment the call to "psplot" if you want contour lines ! of various kinds: solid, dashed, etc ! call psplot ( rhoo, nx, x, nx, y, ny, nlevels, z, xdim, ydim, & ! xs, ys, fileout) CALL atomi (nat, tau, ityp, at, a0, r0, tau1, tau2, xdim, ydim) 20 STOP DEALLOCATE (z) DEALLOCATE (rhoo, x , y ) DEALLOCATE (tau, ityp, z) DEALLOCATE (rhoi, xi, yi) END PROGRAM plotrho ! !----------------------------------------------------------------------- SUBROUTINE cplot (d, imax, jmax, x, xmin, xmax, iub, y, ymin, & ymax, jub, nc, z, xdim, ydim, xs, ys, str, filename) !----------------------------------------------------------------------- ! ! draws a contour plot of d(i,j). PostScript output on unit 1 ! F.Gygi Dec.15 1987 - P. Giannozzi Oct.6 1989 and later ! Algorithm by Paul D. Bourke, Byte magazine, june 1987, p. 143 ! d(0:imax,0:jmax) contains the function to plot ! x(0:imax) workspace ! xmin,xmax determines the range of the variable x ! iub is the number of intervals along the x axis (<=imax) ! y(0:jmax) workspace ! idem for ymin,ymax and jub (<=jmax) ! nc is the number of levels wanted ( <=ncmax), ! z(0:nc) are the levels, ! xdim and ydim are the physical dimensions of the figure in cm ! xs and ys determine a shift of the origin in cm IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER :: imax, jmax, iub, jub, nc real(DP) :: d (0:imax, 0:jmax), x (0:imax), y (0:jmax), z (0:nc) real(DP) :: xmin, xmax, ymin, ymax, xdim, ydim, xs, ys CHARACTER (len=*) :: filename, str INTEGER, PARAMETER :: ncmax = 19 INTEGER :: i, j, k real(DP) :: gray (0:ncmax), dim real(DP) :: cm = 28.453d0, width=0.5d0, gray0=1.0d0, deltagray=0.7d0 ! cm : number of points per cm ! width: linewidth of the contour plot for PostScript printer OPEN (unit = 1, file = filename, status = 'unknown', form = & 'formatted') IF (nc > ncmax .or. nc < 1) STOP ' nc too big or wrong' IF (iub > imax .or.iub < 1) STOP ' iub too big or wrong' IF (jub > jmax .or.jub < 1) STOP ' jub too big or wrong' IF (xdim < 3.0d0 .or. ydim < 3.0d0) STOP ' really too small!' IF (xdim > 20.0.or. ydim > 30.0) STOP ' really too big!' IF (abs (xs) > 20.or. abs (ys) > 30) STOP ' xs or ys are weird' ! initializations for PostScript output WRITE (1, '(a)') '%! PS-Adobe-1.0' WRITE (1, '("%%BoundingBox:",4f6.1)') xs * cm, ys * cm, (xs + & xdim) * cm, (ys + ydim) * cm WRITE (1, '(a)') '/localdict 100 dict def' WRITE (1, '(a)') 'localdict begin' WRITE (1, '(a)') '/cm {28.453 mul} def' WRITE (1, '(a)') '/title {('//str//')} def' WRITE (1, '(a)') '/Times-Roman findfont 12 scalefont setfont' WRITE (1, '(a)') '% cshow prints a centered string at current position' WRITE (1, '(a)') '/cshow {gsave dup stringwidth pop 2 div neg 0' WRITE (1, '(a)') ' rmoveto show grestore} def' WRITE (1, '(a)') '% x1 y1 x2 y2 p : draws a segment from point 1 to point 2' WRITE (1, '(a)') '/p {0 setgray newpath moveto lineto stroke} def' WRITE (1, '(a)') '% x1 y1 x2 y2 x3 y3 x4 y4 sn :' WRITE (1, '(a)') '% fills the region bounded by points 1 to 4' WRITE (1, '(a)') '% with greyscale n' ! type of gray for shaded areas DO k = 0, nc gray (k) = gray0 - k * deltagray / nc IF (k<10) THEN WRITE (1, '("/s",i1," {",f4.2," setgray newpath ", & & "moveto lineto lineto lineto fill} def")') k, gray (k) WRITE (1, '("/t",i1," {",f4.2," setgray newpath ", & & "moveto lineto lineto fill} def")') k, gray (k) ELSE WRITE (1, '("/u",i1," {",f4.2," setgray newpath ", & & "moveto lineto lineto lineto fill} def")') mod (k, 10) , & &gray (k) WRITE (1, '("/v",i1," {",f4.2," setgray newpath ", & & "moveto lineto lineto fill} def")') mod (k, 10) , gray (k & &) ENDIF ENDDO WRITE (1, '(a)') '%%EndPreamble' WRITE (1, '(a)') 'gsave' WRITE (1, '(1x,f6.2," cm ",f6.2," cm translate")') xs, ys WRITE (1, '(a)') '% Uncomment next line if you want a big picture' WRITE (1, '(a)') '% 1.8 1.8 scale' WRITE (1, '(f7.3," setlinewidth")') width WRITE (1, '(a)') '% Comment the next line to remove the title' WRITE (1, '(1x,f6.2," cm ",f6.2," cm moveto title cshow")') & xdim / 2, ydim + 1.5d0 CALL hatch (0.d0, xdim, 0.d0, ydim) DO i = 0, iub x (i) = xdim * dble (i) / iub ENDDO DO j = 0, jub y (j) = ydim * dble (j) / jub ENDDO CALL conrec (imax, iub, jmax, jub, x, y, d, nc, z) ! draw frame of size xdim by ydim WRITE (1, '(a)') '1 setlinewidth 0 setgray newpath' WRITE (1, '(2f6.1," moveto")') 0.0d0, 0.0d0 WRITE (1, '(2f6.1," lineto")') xdim * cm, 0.0d0 WRITE (1, '(2f6.1," lineto")') xdim * cm, ydim * cm WRITE (1, '(2f6.1," lineto")') 0.0d0, ydim * cm WRITE (1, '(a)') 'closepath stroke' ! write (1,'(a)') 'grestore' ! write (1,'(a)') '%%Trailer' ! write (1,'(a)') 'showpage' ! close(1) RETURN END SUBROUTINE cplot SUBROUTINE conrec (imax, iub, jmax, jub, x, y, d, nc, z) IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER :: imax, iub, jmax, jub, nc real(DP) :: d (0:imax, 0:jmax), x (0:imax), y (0:jmax), z (0:nc) INTEGER, PARAMETER :: ncmax = 19 CHARACTER (len=4) :: triangle (0:ncmax), trapez (0:ncmax) real(DP) :: h (0:4), xh (0:4), yh (0:4) real(DP) :: x1, y1, x2, y2, x3, y3, x4, y4, dx, dy, xx, yy, dmin, dmax real(DP) :: cm = 28.453d0 ! cm : conversion factor from cm to points for PostScript INTEGER :: ish (0:4), im (0:3), jm (0:3), castab (0:2, 0:2, 0:2) INTEGER :: i, j, k, m, m1, m2, m3, npoint, icase, levelin, nolevel data (im (i), i = 0, 3) / 0, 1, 1, 0 / data (jm (i), i = 0, 3) / 0, 0, 1, 1 / data ( ( (castab (i, j, k), k = 0, 2), j = 0, 2), i = 0, 2) & / 0, 0, 8, 0, 2, 5, 7, 6, 9, 0, 3, 4, 1, 3, 1, 4, 3, 0, 9, 6, 7, & 5, 2, 0, 8, 0, 0 / dy = (y (jub) - y (0) ) / (nc + 1) xx = x (iub) + 1.0d0 dx = 0.5d0 WRITE (1, '(a)') '% Start of Color Code' CALL hatch (xx, xx + dx, y (0), y (jub) ) DO k = 0, nc yy = y (jub) - k * dy WRITE (1, '(8f6.1,$)') xx * cm, yy * cm, (xx + dx) * cm, yy * cm, & (xx + dx) * cm, (yy - dy) * cm, xx * cm, (yy - dy) * cm IF (k < 10) THEN WRITE (triangle (k) , '(" t",i1,1x)') k WRITE (trapez (k) , '(" s",i1,1x)') k ELSE WRITE (triangle (k) , '(" v",i1)') mod (k, 10) WRITE (trapez (k) , '(" u",i1)') mod (k, 10) ENDIF WRITE (1, '(a4)') trapez (k) WRITE (1, * ) '0 setgray newpath' WRITE (1, '(2f6.1," moveto")') xx * cm, yy * cm WRITE (1, '(2f6.1," lineto")') (xx + dx) * cm, yy * cm WRITE (1, '(2f6.1," lineto")') (xx + dx) * cm, (yy - dy) & * cm WRITE (1, '(2f6.1," lineto")') xx * cm, (yy - dy) * cm WRITE (1, * ) 'closepath stroke' WRITE (1, '(2f6.1," moveto")') (x (iub) + 2.0d0) * cm, (yy - dy / & 2) * cm IF (k == 0) THEN WRITE (1, '("(z<",f7.5,") show")') z (0) ELSEIF (k==nc) THEN WRITE (1, '("(z>",f7.5,") show")') z (nc - 1) ELSE WRITE (1, '("(",f7.5,"<z<",f7.5,") show")') z (k - 1) , & z (k) ENDIF ENDDO WRITE (1, '(a)') '% End of Color Code' DO k = 1, nc - 1 IF (z (k) <=z (k - 1) ) STOP 'zk order' ENDDO ! scan the array, top down, left to right, to paint shaded areas DO j = jub - 1, 0, - 1 DO i = 0, iub - 1 ! find lowest and highest vertex dmin = min (d (i, j), d (i, j + 1), d (i + 1, j), d (i + 1, j + 1) & ) dmax = max (d (i, j), d (i, j + 1), d (i + 1, j), d (i + 1, j + 1) & ) ! search for levels in this box nolevel = 0 DO k = 0, nc - 1 IF (z (k) < dmin) nolevel = k + 1 IF (z (k) >= dmin .and. z (k) <= dmax) THEN levelin = k GOTO 10 ENDIF ENDDO ! no level in this box: paint the whole box and pass to another box WRITE (1, '(8f6.1,a4)') x (i) * cm, y (j) * cm, x (i + 1) * cm, y & (j) * cm, x (i + 1) * cm, y (j + 1) * cm, x (i) * cm, y (j + 1) & * cm, trapez (nolevel) GOTO 100 ! there is at least a level in this box: paint the whole box 10 CONTINUE WRITE (1, '(8f6.1,a4)') x (i) * cm, y (j) * cm, x (i + 1) * cm, y & (j) * cm, x (i + 1) * cm, y (j + 1) * cm, x (i) * cm, y (j + 1) & * cm, trapez (levelin) DO k = levelin, nc - 1 ! if no more levels in this box, move to another box IF (z (k) >dmax) GOTO 100 ! find contour of zero levels in this box DO m = 1, 4 h (m) = d (i + im (m - 1), j + jm (m - 1) ) - z (k) xh (m) = x (i + im (m - 1) ) yh (m) = y (j + jm (m - 1) ) ENDDO h (0) = (h (1) + h (2) + h (3) + h (4) ) / 4 xh (0) = (x (i) + x (i + 1) ) / 2 yh (0) = (y (j) + y (j + 1) ) / 2 DO m = 0, 4 IF (h (m) >0) THEN ish (m) = 2 ELSEIF (h (m) <0) THEN ish (m) = 0 ELSE ish (m) = 1 ENDIF ENDDO ! scan each triangle in the box to paint shaded areas DO m = 1, 4 m1 = m m2 = 0 m3 = mod (m, 4) + 1 npoint = 0 icase = castab (ish (m1), ish (m2), ish (m3) ) IF (icase == 0) THEN IF (ish (m1) ==2) THEN ! paint this triangle if positive x1 = xh (m1) y1 = yh (m1) x2 = xh (m2) y2 = yh (m2) x3 = xh (m3) y3 = yh (m3) npoint = 3 ENDIF ELSEIF (icase == 1) THEN ! line between vertices m1 and m2 x1 = xh (m1) y1 = yh (m1) x2 = xh (m2) y2 = yh (m2) IF (ish (m3) ==2) THEN x3 = xh (m3) y3 = yh (m3) npoint = 3 ENDIF ELSEIF (icase == 2) THEN ! line between vertices m2 and m3 x1 = xh (m2) y1 = yh (m2) x2 = xh (m3) y2 = yh (m3) IF (ish (m1) == 2) THEN x3 = xh (m1) y3 = yh (m1) npoint = 3 ENDIF ELSEIF (icase == 3) THEN ! line between vertices m3 and m1 x1 = xh (m3) y1 = yh (m3) x2 = xh (m1) y2 = yh (m1) IF (ish (m2) == 2) THEN x3 = xh (m2) y3 = yh (m2) npoint = 3 ENDIF ELSEIF (icase == 4) THEN ! line between vertex m1 and side m2-m3 x1 = xh (m1) y1 = yh (m1) x2 = (h (m3) * xh (m2) - h (m2) * xh (m3) ) / (h (m3) - h (m2) & ) y2 = (h (m3) * yh (m2) - h (m2) * yh (m3) ) / (h (m3) - h (m2) & ) IF (ish (m3) == 2) THEN x3 = xh (m3) y3 = yh (m3) ELSE x3 = xh (m2) y3 = yh (m2) ENDIF npoint = 3 ELSEIF (icase == 5) THEN ! line between vertex m2 and side m3-m1 x1 = xh (m2) y1 = yh (m2) x2 = (h (m1) * xh (m3) - h (m3) * xh (m1) ) / (h (m1) - h (m3) & ) y2 = (h (m1) * yh (m3) - h (m3) * yh (m1) ) / (h (m1) - h (m3) & ) IF (ish (m1) == 2) THEN x3 = xh (m1) y3 = yh (m1) ELSE x3 = xh (m3) y3 = yh (m3) ENDIF npoint = 3 ELSEIF (icase == 6) THEN ! line between vertex m3 and line m1-m2 x1 = xh (m3) y1 = yh (m3) x2 = (h (m2) * xh (m1) - h (m1) * xh (m2) ) / (h (m2) - h (m1) & ) y2 = (h (m2) * yh (m1) - h (m1) * yh (m2) ) / (h (m2) - h (m1) & ) IF (ish (m2) == 2) THEN x3 = xh (m2) y3 = yh (m2) ELSE x3 = xh (m1) y3 = yh (m1) ENDIF npoint = 3 ELSEIF (icase == 7) THEN ! line between sides m1-m2 and m2-m3 x1 = (h (m2) * xh (m1) - h (m1) * xh (m2) ) / (h (m2) - h (m1) & ) y1 = (h (m2) * yh (m1) - h (m1) * yh (m2) ) / (h (m2) - h (m1) & ) x2 = (h (m3) * xh (m2) - h (m2) * xh (m3) ) / (h (m3) - h (m2) & ) y2 = (h (m3) * yh (m2) - h (m2) * yh (m3) ) / (h (m3) - h (m2) & ) IF (ish (m2) == 2) THEN x3 = xh (m2) y3 = yh (m2) npoint = 3 ELSE x3 = xh (m3) y3 = yh (m3) x4 = xh (m1) y4 = yh (m1) npoint = 4 ENDIF ELSEIF (icase == 8) THEN ! line between sides m2-m3 and m3-m1 x1 = (h (m3) * xh (m2) - h (m2) * xh (m3) ) / (h (m3) - h (m2) & ) y1 = (h (m3) * yh (m2) - h (m2) * yh (m3) ) / (h (m3) - h (m2) & ) x2 = (h (m1) * xh (m3) - h (m3) * xh (m1) ) / (h (m1) - h (m3) & ) y2 = (h (m1) * yh (m3) - h (m3) * yh (m1) ) / (h (m1) - h (m3) & ) IF (ish (m3) == 2) THEN x3 = xh (m3) y3 = yh (m3) npoint = 3 ELSE x3 = xh (m1) y3 = yh (m1) x4 = xh (m2) y4 = yh (m2) npoint = 4 ENDIF ELSEIF (icase == 9) THEN ! line between sides m3-m1 and m1-m2 x1 = (h (m1) * xh (m3) - h (m3) * xh (m1) ) / (h (m1) - h (m3) & ) y1 = (h (m1) * yh (m3) - h (m3) * yh (m1) ) / (h (m1) - h (m3) & ) x2 = (h (m2) * xh (m1) - h (m1) * xh (m2) ) / (h (m2) - h (m1) & ) y2 = (h (m2) * yh (m1) - h (m1) * yh (m2) ) / (h (m2) - h (m1) & ) IF (ish (m1) == 2) THEN x3 = xh (m1) y3 = yh (m1) npoint = 3 ELSE x3 = xh (m2) y3 = yh (m2) x4 = xh (m3) y4 = yh (m3) npoint = 4 ENDIF ENDIF IF (npoint == 3) THEN WRITE (1, '(6f6.1,a4)') x1 * cm, y1 * cm, x2 * cm, y2 * cm, x3 & * cm, y3 * cm, triangle (k + 1) ELSEIF (npoint == 4) THEN WRITE (1, '(8f6.1,a4)') x1 * cm, y1 * cm, x2 * cm, y2 * cm, x3 & * cm, y3 * cm, x4 * cm, y4 * cm, trapez (k + 1) ENDIF IF (icase/=0) WRITE (1, '(4f6.1," p")') x1 * cm, y1 * cm, x2 * & cm, y2 * cm ENDDO ENDDO 100 CONTINUE ENDDO ENDDO ! RETURN END SUBROUTINE conrec ! !----------------------------------------------------------------------- SUBROUTINE atomi (nat, tau, ityp, at, a0, r0, tau1, tau2, xdim, ydim) !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER :: nat, ityp (nat) real(DP) :: tau (3, nat), at (3, 3), r0 (3), tau1 (3), tau2 (3), a0 real(DP) :: xdim, ydim INTEGER :: n1, n2, n3, i, n real(DP) :: r (3), ri (3), tau1n, tau2n, delta0, r1, r2, r3 real(DP) :: delta = 1.0d0, cm = 28.453d0 ! delta0 = delta / a0 tau1n = sqrt (tau1 (1) **2 + tau1 (2) **2 + tau1 (3) **2) tau2n = sqrt (tau2 (1) **2 + tau2 (2) **2 + tau2 (3) **2) ! clip the rectangle WRITE (1, '(a)') 'gsave newpath' WRITE (1, '(2f6.1," moveto")') 0.0d0, 0.0d0 WRITE (1, '(2f6.1," lineto")') xdim * cm, 0.0d0 WRITE (1, '(2f6.1," lineto")') xdim * cm, ydim * cm WRITE (1, '(2f6.1," lineto")') 0.0d0, ydim * cm WRITE (1, '(a)') 'closepath clip stroke' ! really bad algorithm to generate (hopefully) all atoms ! that are inside the contour plane or at |z| < delta (a.u.) DO n1 = - 3, + 3 DO n2 = - 3, + 3 DO n3 = - 3, + 3 DO i = 1, 3 r (i) = n1 * at (i, 1) + n2 * at (i, 2) + n3 * at (i, 3) ENDDO DO n = 1, nat DO i = 1, 3 ri (i) = tau (i, n) + r (i) - r0 (i) ENDDO ! this is the component along direction 1 ... r1 = (ri (1) * tau1 (1) + ri (2) * tau1 (2) + ri (3) * tau1 (3))& / tau1n IF (r1 > - delta0 .and. r1 < tau1n + delta0) THEN ! and this is along direction 2 ... r2 = (ri (1) * tau2(1) + ri (2) * tau2(2) + ri (3) * tau2(3))& / tau2n IF (r2 > - delta0 .and. r2 < tau2n + delta0) THEN DO i = 1, 3 ri (i) = ri (i) - r1 * tau1 (i) / tau1n & - r2 * tau2 (i) / tau2n ENDDO r3 = sqrt (ri (1) **2 + ri (2) **2 + ri (3) **2) ! and this is along the direction orthogonal to plane IF (abs (r3) < delta0) THEN WRITE (1, '(3f6.1," 0 360 arc gsave ",f4.2, & & " setgray fill grestore stroke")') & r1 / tau1n * xdim * cm, r2 / tau2n * ydim * cm,& delta0 / tau1n * xdim * cm, abs (r3) / delta0 ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO ! WRITE (1, '(a)') 'grestore' WRITE (1, '(a)') 'grestore' WRITE (1, '(a)') '%%Trailer' WRITE (1, '(a)') 'showpage' CLOSE (1) RETURN END SUBROUTINE atomi SUBROUTINE hatch (x1, x2, y1, y2) IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) real(DP) :: x1, x2, y1, y2 real(DP) :: cm = 28.453d0, delta = 0.2d0, dim INTEGER :: nhach, n WRITE (1, '(a)') '% Beginning of hatching' WRITE (1, '(a)') 'gsave newpath' WRITE (1, '(2f6.1," moveto")') x1 * cm, y1 * cm WRITE (1, '(2f6.1," lineto")') x2 * cm, y1 * cm WRITE (1, '(2f6.1," lineto")') x2 * cm, y2 * cm WRITE (1, '(2f6.1," lineto")') x1 * cm, y2 * cm WRITE (1, '(a)') 'closepath clip' dim = max (x2 - x1, y2 - y1) nhach = dim / delta ! delta=dim/nhach DO n = 1, nhach WRITE (1, '(4f6.1," p")') (x1 + (n - 1) * delta) * cm, y1 * cm, & (x1 + dim) * cm, (y1 + dim - (n - 1) * delta) * cm WRITE (1, '(4f6.1," p")') x1 * cm, (y1 + (n - 1) * delta) & * cm, (x1 + dim - (n - 1) * delta) * cm, (y1 + dim) * cm ENDDO WRITE (1, '(a)') 'grestore' WRITE (1, '(a)') '% End of hatching' RETURN END SUBROUTINE hatch ! !----------------------------------------------------------------------- SUBROUTINE psplot ( d, imax, x, iub, y, jub, nlevels, z, & xdim, ydim, xs, ys, filename) !----------------------------------------------------------------------- ! ! draws a contour plot of d(i,j). PostScript output on unit 1 ! ! d(0:imax,0:jub) contains the function to plot on a uniform 2d grid ! x(0:iub) workspace ! iub is the number of intervals along the x axis (<=imax) ! y(0:jub) workspace ! jub is the number of intervals along the y axis ! nlevels is the number of levels wanted ! z(0:nlevels) are the levels ! xdim and ydim are the physical dimensions of the figure in cm ! xs and ys determine a shift of the origin in cm IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER, PARAMETER :: ncontourmax=100, npointmax=500, nmaxtypes=8 INTEGER :: imax, iub, jub, nlevels real(DP) :: d(0:imax,0:jub), x(0:imax), y(0:jub), z(0:nlevels), & xdim, ydim, xs, ys ! real(DP) :: line(2,npointmax,ncontourmax), work(2,npointmax), & segments(2,2,npointmax*ncontourmax), yy, dy INTEGER :: npoints(ncontourmax), ncontours, i, j, k, l, npts real(DP), PARAMETER :: cm = 28.453d0 LOGICAL :: lwork(npointmax*ncontourmax) CHARACTER (len=256) :: filename, linetype(nmaxtypes) data linetype & / '[1 2 ] 0 setdash 0.5 setlinewidth', & '[1 2 4 2] 0 setdash', & '[4 2 ] 0 setdash', & '[ ] 0 setdash', & '[1 2 ] 0 setdash 1.0 setlinewidth', & '[1 2 4 2] 0 setdash', & '[4 2 ] 0 setdash', & '[ ] 0 setdash'/ OPEN (unit = 1, file = filename, status = 'unknown', form = & 'formatted') IF(nlevels < 1) STOP ' nunmber of levels is wrong' IF(iub > imax .or. iub < 1) STOP ' iub too big or wrong' IF(jub < 1) STOP ' jub or wrong' IF(xdim < 3.0d0 .or. ydim < 3.0d0 ) STOP ' really too small!' IF(xdim > 20.0.or. ydim > 30.0) STOP ' really too big!' IF(abs(xs) > 20 .or. abs(ys) > 30) STOP ' xs or ys are weird' ! initializations for PostScript output WRITE (1,'(a)') '%! PS-Adobe-1.0' WRITE (1,'(''%%BoundingBox:'',4f6.1)') & xs*cm, ys*cm,(xs+xdim)*cm,(ys+ydim)*cm WRITE (1,'(a)') '/localdict 100 dict def' WRITE (1,'(a)') 'localdict begin' WRITE (1,'(a)') '/cm {28.453 mul} def' WRITE (1,'(a)') '/title {('//filename//')} def' WRITE (1,'(a)') '/Times-Roman findfont 12 scalefont setfont' WRITE (1,'(a)') '% cshow prints a centered string at current position' WRITE (1,'(a)') '/cshow {gsave dup stringwidth pop 2 div neg 0' WRITE (1,'(a)') ' rmoveto show grestore} def' WRITE (1,'(a)') '% x1 y1 nm x2 y2 ...xn yn n-1 lrs :' WRITE (1,'(a)') '% draws a line from point 1 to point n' WRITE (1,'(a)') '/nm {newpath moveto} def' WRITE (1,'(a)') '/lrs {{lineto} repeat stroke} def' WRITE (1,'(a)') '%%EndPreamble' WRITE (1,'(a)') 'gsave' WRITE (1,'(a)') '0 setgray' WRITE (1,'(1x,f6.2,'' cm '',f6.2,'' cm translate'')') xs, ys WRITE (1,'(1x,f6.2,'' cm '',f6.2,'' cm moveto title cshow'')') & xdim/2, ydim+1.5d0 ! DO i=0,iub x(i)=xdim*dble(i)/iub ENDDO DO j=0,jub y(j)=ydim*dble(j)/jub ENDDO DO k=1,nlevels-1 IF (z(k)<=z(k-1)) STOP ' order of levels' ENDDO yy=y(jub)*0.75d0 dy=(y(jub)-y(0))/(nlevels-1)/2 DO k=0,nlevels-1 ! type of line WRITE (1,'(''%% This is level '',i3)') k WRITE (1,'(50a)') linetype(mod(k,nmaxtypes)+1) WRITE (1,*) ' newpath ' WRITE (1,'(2f6.1,'' moveto'')') (x(iub)+1.0d0)*cm, yy*cm WRITE (1,'(2f6.1,'' lineto'')') (x(iub)+1.8d0)*cm, yy*cm WRITE (1,*) 'closepath stroke' WRITE (1,'(2f6.1,'' moveto'')') (x(iub)+2.0d0)*cm, yy*cm WRITE (1,'(''( z = '',f6.3,'') show'')') z(k) yy=yy-dy CALL findcontours ( d, imax, iub, jub, x, y, z(k), segments, work, & lwork, npointmax, ncontourmax, line, npoints, ncontours) DO l=1,ncontours WRITE (1,'(12f6.1)') & ( line(1,npts,l)*cm, line(2,npts,l)*cm, npts=npoints(l),1,-1 ) WRITE (1,'(''nm'',i4,'' lrs'')') npoints(l)-1 ENDDO ENDDO ! draw frame of size xdim by ydim WRITE (1,*) '[] 0 setdash 1 setlinewidth 0 setgray newpath' WRITE (1,'(2f6.1,'' moveto'')') 0.0d0,0.0d0 WRITE (1,'(2f6.1,'' lineto'')') xdim*cm,0.0d0 WRITE (1,'(2f6.1,'' lineto'')') xdim*cm,ydim*cm WRITE (1,'(2f6.1,'' lineto'')') 0.0d0,ydim*cm WRITE (1,*) 'closepath stroke' ! write (1,*) 'grestore' ! write (1,'(a)') '%%Trailer' ! write (1,*) 'showpage' RETURN END SUBROUTINE psplot ! SUBROUTINE findcontours (d, imax, iub, jub, x, y, z, segments, work, & done, npointmax, ncontourmax, line, npoints, ncontours) ! IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER :: imax, iub, jub, npointmax, ncontourmax, ncontours, & npoints(ncontourmax) real(DP) :: d(0:imax,0:jub),x(0:iub), y(0:jub), z, & line(2,npointmax,ncontourmax), segments(2,2,npointmax*ncontourmax) LOGICAL :: done(npointmax*ncontourmax) ! INTEGER :: i, j, m, i0, j0, m0, nsegments, nsegmax, nseg, nnext, npts, & isame, inext, nleft real(DP) :: epsx, epsy, work(2,npointmax) LOGICAL :: found, start_from_boundary, hit_boundary, loop_closed nsegmax=npointmax*ncontourmax nsegments=0 DO j0=jub-1,0,-1 DO i0=0,iub-1 DO m0=1,4 ! scan boxes and triangles until a contour is found CALL triangle ( d, x, y, z, imax, iub, jub, & i0, j0, m0, found, segments(1,1,nsegments+1) ) IF (found) THEN nsegments=nsegments+1 IF (nsegments > nsegmax) THEN PRINT '('' Erore: nsegments, nsegmax '',2i5)', & nsegments, nsegmax STOP ENDIF ENDIF ENDDO ENDDO ENDDO DO nseg=1,nsegments done(nseg)=.false. ENDDO ncontours=0 nleft=nsegments epsx = 0.0001d0 * (x(1)-x(0)) epsy = 0.0001d0 * (y(1)-y(0)) DO nseg=1,nsegments IF (.not.done(nseg)) THEN ncontours=ncontours+1 IF (ncontours>ncontourmax) THEN PRINT '('' too many ('',i4,'') contours'')', ncontours STOP ENDIF line(1,1,ncontours)=segments(1,1,nseg) line(2,1,ncontours)=segments(2,1,nseg) start_from_boundary = & abs(segments(1,1,nseg)-x( 0)) < epsx .or. & abs(segments(1,1,nseg)-x(iub)) < epsx .or. & abs(segments(2,1,nseg)-y( 0)) < epsy .or. & abs(segments(2,1,nseg)-y(jub)) < epsy ! true if the starting point lies on a boundary nnext=nseg inext=2 npts=1 10 npts=npts+1 IF (npts > npointmax) THEN PRINT '('' too many ('',i5,'') points'')', npts STOP ENDIF line(1,npts,ncontours)=segments(1,inext,nnext) line(2,npts,ncontours)=segments(2,inext,nnext) nleft=nleft-1 done(nnext)=.true. ! check if the loop is closed loop_closed = ( abs(segments(1,inext,nnext) - & line(1,1,ncontours)) < epsx ) .and. & ( abs(segments(2,inext,nnext) - & line(2,1,ncontours)) < epsy ) ! check if the border has been reached hit_boundary = & abs(segments(1,inext,nnext)-x( 0)) < epsx .or. & abs(segments(1,inext,nnext)-x(iub)) < epsx .or. & abs(segments(2,inext,nnext)-y( 0)) < epsy .or. & abs(segments(2,inext,nnext)-y(jub)) < epsy IF (nleft == 0) THEN IF (.not.hit_boundary .and. .not.loop_closed) & PRINT '('' Warning: no more points'')' npoints(ncontours)=npts RETURN ENDIF IF (start_from_boundary) THEN IF (hit_boundary) THEN npoints(ncontours)=npts GOTO 20 ELSEIF (loop_closed) THEN PRINT '('' Warning: loop closed on the boundary ?'')' GOTO 20 ENDIF ELSE IF (loop_closed) THEN npoints(ncontours)=npts GOTO 20 ELSEIF (hit_boundary) THEN DO nnext=1,npts work(1,nnext)=line(1,npts-nnext+1,ncontours) work(2,nnext)=line(2,npts-nnext+1,ncontours) ENDDO DO nnext=1,npts line(1,nnext,ncontours)=work(1,nnext) line(2,nnext,ncontours)=work(2,nnext) ENDDO start_from_boundary=.true. ENDIF ENDIF DO nnext=1,nsegments IF (.not.done(nnext)) THEN DO isame=1,2 found= ( abs(segments(1,isame,nnext) - & line(1,npts,ncontours)) < epsx ) .and. & ( abs(segments(2,isame,nnext) - & line(2,npts,ncontours)) < epsy ) IF (found) THEN inext=mod(isame,2)+1 !inext =1 if isame=2 and vice versa GOTO 10 ENDIF ENDDO ENDIF ENDDO npoints(ncontours)=npts-1 PRINT '('' Warning: next point not found'',i5)', nseg ENDIF 20 CONTINUE ENDDO RETURN END SUBROUTINE findcontours ! SUBROUTINE triangle ( d, x, y, z, imax, iub, jub, i0, j0, m0, found, line ) ! IMPLICIT NONE INTEGER, PARAMETER :: DP = selected_real_kind(14,200) INTEGER :: imax, iub, jub, i0, j0, m0 real(DP) :: d(0:imax,0:jub), x(0:iub), y(0:jub), z, line(2,2) LOGICAL :: found ! real(DP) :: h(0:4), xh(0:4), yh(0:4), dmin, dmax INTEGER :: m, icase, m1, m2, m3, i, j, k, ish(0:4), castab(-1:1,-1:1,-1:1) data (((castab(i,j,k),k=-1,1),j=-1,1),i=-1,1) / & 0 , 0 , 8 , 0 , 2 , 5 , 7 , 6 , 9 , 0 , 3 , 4 , & 1 , 3 , 1 , 4 , 3 , 0 , 9 , 6 , 7 , 5 , 2 , 0 , 8 , 0 , 0 / ! values of the function in the four corners... h (1) = d(i0 ,j0 ) - z h (2) = d(i0+1,j0 ) - z h (3) = d(i0+1,j0+1) - z h (4) = d(i0 ,j0+1) - z dmin=min(h(1),h(2),h(3),h(4)) dmax=max(h(1),h(2),h(3),h(4)) ! return if no contour is found in this box IF ( dmin > 0.0d0 .or. dmax < 0.0d0 ) THEN found=.false. RETURN ENDIF ! xy coordinates of the four corners... xh(1) = x(i0 ) xh(2) = x(i0+1) xh(3) = x(i0+1) xh(4) = x(i0 ) yh(1) = y(j0 ) yh(2) = y(j0 ) yh(3) = y(j0+1) yh(4) = y(j0+1) ! and this is the center of the box h(0) =(h(1)+h(2)+h(3)+h(4))/4.0d0 xh(0)=(x(i0)+x(i0+1))/2.0d0 yh(0)=(y(j0)+y(j0+1))/2.0d0 ! ish=+1 if h>0 ; ish=0 if h=0 ; ish=-1 if h<0 DO m=0,4 ish(m) = nint( sign(1.d0, h(m)) ) IF (h(m) == 0.0d0) ish(m) = 0 ENDDO ! starting triangle m1=m0 m2=0 m3=mod(m0,4)+1 icase=castab(ish(m1),ish(m2),ish(m3)) IF (icase == 0) THEN ! no contour in this triangle found=.false. RETURN ! NOTA BENE: if there is a point on the boundary, ! this is always chosen to be the first point ELSEIF (icase == 1) THEN ! line between vertices m1 and m2 line(1,1)=xh(m1) line(2,1)=yh(m1) line(1,2)=xh(m2) line(2,2)=yh(m2) ELSEIF (icase == 2) THEN ! line between vertices m2 and m3 line(1,1)=xh(m3) line(2,1)=yh(m3) line(1,2)=xh(m2) line(2,2)=yh(m2) ELSEIF (icase == 3) THEN ! line between vertices m3 and m1 line(1,1)=xh(m3) line(2,1)=yh(m3) line(1,2)=xh(m1) line(2,2)=yh(m1) ELSEIF (icase == 4) THEN ! line between vertex m1 and side m2-m3 line(1,1)=xh(m1) line(2,1)=yh(m1) line(1,2)=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) line(2,2)=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) ELSEIF (icase == 5) THEN ! line between vertex m2 and side m3-m1 line(1,1)=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) line(2,1)=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) line(1,2)=xh(m2) line(2,2)=yh(m2) ELSEIF (icase == 6) THEN ! line between vertex m3 and line m1-m2 line(1,1)=xh(m3) line(2,1)=yh(m3) line(1,2)=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) line(2,2)=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) ELSEIF (icase == 7) THEN ! line between sides m1-m2 and m2-m3 line(1,1)=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) line(2,1)=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) line(1,2)=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) line(2,2)=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) ELSEIF (icase == 8) THEN ! line between sides m2-m3 and m3-m1 line(1,1)=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) line(2,1)=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) line(1,2)=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) line(2,2)=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) ELSEIF (icase == 9) THEN ! line between sides m3-m1 and m1-m2 line(1,1)=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) line(2,1)=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) line(1,2)=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) line(2,2)=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) ENDIF found=.true. RETURN END SUBROUTINE triangle
_______________________________________________ Pw_forum mailing list Pw_forum@pwscf.org http://pwscf.org/mailman/listinfo/pw_forum
