scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Chnages proposed by Leonard Talirz to allow multiple plotting 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@13223 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2017-01-05 13:55:02 +00:00
parent ae5c7bad26
commit 7570c21f63
8 changed files with 573 additions and 388 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
PP/Doc/INPUT_PP.def
View File

@ -6,8 +6,8 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
@b {Purpose of pp.x:} data analysis and plotting.
The code performs two steps:
(1) reads the output file produced by pw.x, extract and calculate
the desired quantity (rho, V, ...)
(1) reads the output produced by pw.x, extracts and calculates
the desired quantity/quantities (rho, V, ...)
(2) writes the desired quantity to file in a suitable format for
various types of plotting and various plotting programs
@ -71,7 +71,7 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
2 = local ionic potential V_bare
3 = local density of states at e_fermi
3 = local density of states at specific energy or grid of energies
(number of states per volume, in bohr^3,
per energy unit, in Ry)
@ -91,7 +91,7 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
8 = electron localization function (ELF)
9 = charge density minus superposition of atomic densities
9 = charge density minus superposition of atomic densities
10 = integrated local density of states (ILDOS)
from @ref emin to @ref emax (emin, emax in eV)
@ -107,15 +107,15 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
can be performed for PAW calculations only
requires a very dense real-space grid!
18 = The exchange and correlation magnetic field in
18 = The exchange and correlation magnetic field in
the noncollinear case
19 = Reduced density gradient
19 = Reduced density gradient
(J. Chem. Theory Comput. 7, 625 (2011))
Set the isosurface between 0.3 and 0.6 to plot the
non-covalent interactions (see also plot_num = 20)
20 = Product of the electron density (charge) and the second
20 = Product of the electron density (charge) and the second
eigenvalue of the electron-density Hessian matrix;
used to colorize the RDG plot (plot_num = 19)
@ -156,6 +156,51 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
}
}
}
elsewhen -test "plot_num=3" {
label {
Options for LDOS (plot_num=3):
LDOS is plotted on grid [emin, emax] with spacing delta_e.
}
var emin -type REAL {
default e_fermi
info {
lower boundary of energy grid (in eV).
Defaults to Fermi energy.
}
}
var emax -type REAL {
status OPTIONAL
info {
upper boundary of energy grid (in eV).
If not specified, LDOS is computed just for energy @ref emin
}
}
var delta_e -type REAL {
default 0.1
status OPTIONAL
info {
spacing of energy grid (in eV).
i.e. compute ILDOS from @ref emin to @ref emax
}
}
var degauss_ldos -type REAL {
default {degauss (converted to eV)}
status OPTIONAL
info {
broadening of energy levels for LDOS (in eV).
Defaults to broadening specified for electronic smearing
of the calculation.
}
}
}
elsewhen -test "plot_num=5" {
label {
@ -174,18 +219,23 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
label {
Options for |psi|^2 (plot_num=7):
}
var kpoint -type INTEGER {
dimension kpoint -start 1 -end 2 -type INTEGER {
info {
Unpolarized and noncollinear case: k-point to be plotted
LSDA: k-point and spin polarization to be plotted
(spin-up and spin-down correspond to different k-points!)
When both kpoint(1) and kpoint(2) are specified, all
kpoints in the range [kpoint(1),kpoint(2)] are plotted.
}
}
var kband -type INTEGER {
dimension kband -start 1 -end 2 -type INTEGER {
info {
band to be plotted
band to be plotted .
When both kband(1) and kband(2) are specified, all
bands in the range [kband(1),kband(2)] are plotted.
}
}
@ -195,7 +245,9 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
}
}
var spin_component -type INTEGER {
dimension spin_component -start 1 -end 2 -type INTEGER {
default 0
status OPTIONAL
info {
Noncollinear case only:
plot the contribution of the given state to the charge
@ -207,6 +259,10 @@ input_description -distribution {Quantum Espresso} -package PWscf -program pp.x
3 = z.
Ignored in unpolarized or LSDA case
When both spin_component(1) and spin_component(2) are specified, all
components in the range [spin_component(1),spin_component(2)] are
plotted.
}
}
}

2
PP/src/Makefile
View File

@ -14,7 +14,7 @@ add_shift_lc.o \
add_shift_us.o \
atomic_wfc_nc_proj.o \
cft.o \
chdens.o \
chdens_module.o \
chdens_bspline.o \
compute_ppsi.o \
compute_sigma_avg.o \

1
PP/src/chdens_bspline.f90
View File

@ -255,6 +255,7 @@ SUBROUTINE plot_3d_bspline (alat, at, nat, tau, atm, ityp, rhor, &
USE kinds, ONLY : dp
USE io_global, ONLY : stdout, ionode
USE fft_base, ONLY : dfftp
USE chdens_module, ONLY : write_openmol_file
!---------------------------------------------------------------------
implicit none
integer, intent(in) :: nx, ny, nz, nat, ityp(nat), output_format, ounit

550
PP/src/chdens.f90 → PP/src/chdens_module.f90
View File

@ -7,7 +7,10 @@
!
!
!-----------------------------------------------------------------------
SUBROUTINE chdens (filplot,plot_num)
MODULE chdens_module
CONTAINS
!-----------------------------------------------------------------------
SUBROUTINE chdens (plot_files,plot_num)
!-----------------------------------------------------------------------
! Writes the charge density (or potential, or polarisation)
! into a file format suitable for plotting
@ -40,7 +43,8 @@ SUBROUTINE chdens (filplot,plot_num)
USE wavefunctions_module, ONLY: psic
IMPLICIT NONE
CHARACTER (len=256), INTENT(in) :: filplot
!CHARACTER (len=256), INTENT(in) :: filplot
CHARACTER (len=256), DIMENSION(:), ALLOCATABLE, INTENT(in) :: plot_files
!
! If plot_num=-1 the dimensions and structural data are read from the charge
! or potential file, otherwise it uses the data already read from
@ -52,14 +56,14 @@ SUBROUTINE chdens (filplot,plot_num)
! maximum number of files with charge
INTEGER :: ounit, iflag, ios, ipol, nfile, ifile, nx, ny, nz, &
na, i, output_format, idum, direction
na, i, output_format, idum, direction, iplot
real(DP) :: e1(3), e2(3), e3(3), x0 (3), radius, m1, m2, m3, &
weight (nfilemax), isovalue,heightmin,heightmax
real(DP), ALLOCATABLE :: aux(:)
CHARACTER (len=256) :: fileout
CHARACTER (len=256) :: fileout, fileout_tmp
CHARACTER (len=13), DIMENSION(0:7) :: formatname = &
(/ 'gnuplot ', &
'contour.x ', &
@ -98,7 +102,7 @@ SUBROUTINE chdens (filplot,plot_num)
! set the DEFAULT values
!
nfile = 1
filepp(1) = filplot
filepp(1) = plot_files(1)
weight(1) = 1.0d0
iflag = 0
radius = 1.0d0
@ -138,6 +142,7 @@ SUBROUTINE chdens (filplot,plot_num)
RETURN
ENDIF
CALL mp_bcast( filepp, ionode_id, world_comm )
CALL mp_bcast( weight, ionode_id, world_comm )
CALL mp_bcast( iflag, ionode_id, world_comm )
@ -166,7 +171,12 @@ SUBROUTINE chdens (filplot,plot_num)
! check for number of files
!
IF (nfile < 1 .or. nfile > nfilemax) &
CALL errore ('chdens ', 'nfile is wrong ', 1)
CALL errore ('chdens ', 'nfile < 1 or too large', 1)
IF (nfile > 1 .AND. SIZE(plot_files) > 1) THEN
CALL errore ('chdens ', &
"can't mix nfile > 1 with specifying multiple plots", 1)
ENDIF
! check for iflag
@ -279,270 +289,282 @@ SUBROUTINE chdens (filplot,plot_num)
CALL fft_type_allocate ( dffts, at, bg, gcutms, intra_bgrp_comm)
ENDIF
ALLOCATE (rhor(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x))
ALLOCATE (rhos(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x))
ALLOCATE (taus( 3 , nat))
ALLOCATE (ityps( nat))
!
rhor (:) = 0.0_DP
!
! Read files, verify consistency
! Note that only rho is read; all other quantities are discarded
!
DO ifile = 1, nfile
!
CALL plot_io (filepp (ifile), title, nr1sxa, nr2sxa, nr3sxa, &
nr1sa, nr2sa, nr3sa, nats, ntyps, ibravs, celldms, ats, gcutmsa, &
duals, ecuts, idum, atms, ityps, zvs, taus, rhos, - 1)
fileout_tmp = fileout
! Plotting all files in plot_files
DO iplot=1, SIZE(plot_files)
IF (ifile==1.and.plot_num==-1) THEN
atm=atms
ityp=ityps
zv=zvs
tau=taus
ENDIF
!
IF (nats>nat) CALL errore ('chdens', 'wrong file order? ', 1)
IF (dfftp%nr1x/=nr1sxa.or.dfftp%nr2x/=nr2sxa) CALL &
errore ('chdens', 'incompatible nr1x or nr2x', 1)
IF (dfftp%nr1/=nr1sa.or.dfftp%nr2/=nr2sa.or.dfftp%nr3/=nr3sa) CALL &
errore ('chdens', 'incompatible nr1 or nr2 or nr3', 1)
IF (ibravs/=ibrav) CALL errore ('chdens', 'incompatible ibrav', 1)
IF (abs(gcutmsa-gcutm)>1.d-8.or.abs(duals-dual)>1.d-8.or.&
abs(ecuts-ecutwfc)>1.d-8) &
CALL errore ('chdens', 'incompatible gcutm or dual or ecut', 1)
IF (ibravs /= 0 ) THEN
DO i = 1, 6
IF (abs( celldm (i)-celldms (i) ) > 1.0d-7 ) &
CALL errore ('chdens', 'incompatible celldm', 1)
ENDDO
ENDIF
!
rhor (:) = rhor (:) + weight (ifile) * rhos (:)
ENDDO
DEALLOCATE (ityps)
DEALLOCATE (taus)
DEALLOCATE (rhos)
!
! open output file, i.e., "fileout"
!
IF (ionode) THEN
IF (fileout /= ' ') THEN
ounit = 1
OPEN (unit=ounit, file=fileout, form='formatted', status='unknown')
WRITE( stdout, '(/5x,"Writing data to be plotted to file ",a)') &
trim(fileout)
ELSE
ounit = 6
ENDIF
ENDIF
! the isostm subroutine is called only when isostm_flag is true and the
! charge density is related to an STM image (5) or is read from a file
IF ( (isostm_flag) .AND. ( (plot_num == -1) .OR. (plot_num == 5) ) ) THEN
IF ( .NOT. (iflag == 2))&
CALL errore ('chdens', 'isostm should have iflag = 2', 1)
CALL isostm_plot(rhor, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, &
isovalue, heightmin, heightmax, direction)
END IF
IF (SIZE(plot_files) > 1) THEN
filepp(1) = plot_files(iplot)
WRITE(fileout,"(A,A)") TRIM(plot_files(iplot)), TRIM(fileout_tmp)
ENDIF
ALLOCATE (rhor(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x))
ALLOCATE (rhos(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x))
ALLOCATE (taus( 3 , nat))
ALLOCATE (ityps( nat))
!
rhor (:) = 0.0_DP
!
! Read files, verify consistency
! Note that only rho is read; all other quantities are discarded
!
DO ifile = 1, nfile
!
CALL plot_io (filepp (ifile), title, nr1sxa, nr2sxa, nr3sxa, &
nr1sa, nr2sa, nr3sa, nats, ntyps, ibravs, celldms, ats, gcutmsa, &
duals, ecuts, idum, atms, ityps, zvs, taus, rhos, - 1)
!
! At this point we start the calculations, first we normalize the
! vectors defining the plotting region.
! If these vectors have 0 length, replace them with crystal axis
!
m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
IF (abs(m1) < 1.d-6) THEN
e1 (:) = at(:,1)
m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
ENDIF
e1 (:) = e1 (:) / m1
!
m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
IF (abs(m2) < 1.d-6) THEN
e2 (:) = at(:,2)
m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
ENDIF
e2 (:) = e2 (:) / m2
!
m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
IF (abs(m3) < 1.d-6) THEN
e3 (:) = at(:,3)
m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
ENDIF
e3 (:) = e3 (:) / m3
!
! are vectors defining the plotting region aligned along xyz ?
!
fast3d = ( e1(2) == 0.d0 .and. e1(3) == 0.d0) .and. &
( e2(1) == 0.d0 .and. e2(3) == 0.d0) .and. &
( e3(1) == 0.d0 .and. e3(2) == 0.d0)
!
! are crystal axis aligned along xyz ?
!
fast3d = fast3d .and. &
( at(2,1) == 0.d0 .and. at(3,1) == 0.d0) .and. &
( at(1,2) == 0.d0 .and. at(3,2) == 0.d0) .and. &
( at(1,3) == 0.d0 .and. at(2,3) == 0.d0)
fast3d = fast3d .and. (trim(interpolation) == 'fourier')
!
! Initialise FFT for rho(r) => rho(G) conversion if needed
!
IF (.not. ( iflag == 3 .and. ( output_format == 5 .or. &
output_format == 6 .or. &
fast3d ) ) ) THEN
IF (plot_num==-1) THEN
!
gamma_only=.false.
! nproc_pool=1
!
CALL data_structure ( gamma_only )
CALL allocate_fft()
!
! and rebuild G-vectors in reciprocal space
!
CALL ggen ( gamma_only, at, bg )
!
! here we compute the fourier components of the quantity to plot
!
ELSE
!
IF (gamma_only .and. (trim(interpolation) == 'fourier')) THEN
WRITE(stdout,'(/"BEWARE: plot requiring G-space interpolation",&
&" not implemented for Gamma only!",/, &
&"SOLUTION: restart this calculation with", &
&" emtpy namelist &inputpp")')
CALL errore ('chdens','Not implemented, please read above',1)
ENDIF
!
ENDIF
IF (ifile==1.and.plot_num==-1) THEN
atm=atms
ityp=ityps
zv=zvs
tau=taus
ENDIF
!
IF (nats>nat) CALL errore ('chdens', 'wrong file order? ', 1)
IF (dfftp%nr1x/=nr1sxa.or.dfftp%nr2x/=nr2sxa) CALL &
errore ('chdens', 'incompatible nr1x or nr2x', 1)
IF (dfftp%nr1/=nr1sa.or.dfftp%nr2/=nr2sa.or.dfftp%nr3/=nr3sa) CALL &
errore ('chdens', 'incompatible nr1 or nr2 or nr3', 1)
IF (ibravs/=ibrav) CALL errore ('chdens', 'incompatible ibrav', 1)
IF (abs(gcutmsa-gcutm)>1.d-8.or.abs(duals-dual)>1.d-8.or.&
abs(ecuts-ecutwfc)>1.d-8) &
CALL errore ('chdens', 'incompatible gcutm or dual or ecut', 1)
IF (ibravs /= 0 ) THEN
DO i = 1, 6
IF (abs( celldm (i)-celldms (i) ) > 1.0d-7 ) &
CALL errore ('chdens', 'incompatible celldm', 1)
ENDDO
ENDIF
!
rhor (:) = rhor (:) + weight (ifile) * rhos (:)
ENDDO
DEALLOCATE (ityps)
DEALLOCATE (taus)
DEALLOCATE (rhos)
!
! open output file, i.e., "fileout"
!
IF (ionode) THEN
IF (fileout /= ' ') THEN
ounit = 1
OPEN (unit=ounit, file=fileout, form='formatted', status='unknown')
WRITE( stdout, '(/5x,"Writing data to be plotted to file ",a)') &
trim(fileout)
ELSE
ounit = 6
ENDIF
ENDIF
! the isostm subroutine is called only when isostm_flag is true and the
! charge density is related to an STM image (5) or is read from a file
IF ( (isostm_flag) .AND. ( (plot_num == -1) .OR. (plot_num == 5) ) ) THEN
IF ( .NOT. (iflag == 2))&
CALL errore ('chdens', 'isostm should have iflag = 2', 1)
CALL isostm_plot(rhor, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, &
isovalue, heightmin, heightmax, direction)
END IF
!
! At this point we start the calculations, first we normalize the
! vectors defining the plotting region.
! If these vectors have 0 length, replace them with crystal axis
!
m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
IF (abs(m1) < 1.d-6) THEN
e1 (:) = at(:,1)
m1 = sqrt (e1 (1)**2 + e1 (2)**2 + e1 (3)**2)
ENDIF
e1 (:) = e1 (:) / m1
!
m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
IF (abs(m2) < 1.d-6) THEN
e2 (:) = at(:,2)
m2 = sqrt (e2 (1)**2 + e2 (2)**2 + e2 (3)**2)
ENDIF
e2 (:) = e2 (:) / m2
!
m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
IF (abs(m3) < 1.d-6) THEN
e3 (:) = at(:,3)
m3 = sqrt (e3 (1)**2 + e3 (2)**2 + e3 (3)**2)
ENDIF
e3 (:) = e3 (:) / m3
!
! are vectors defining the plotting region aligned along xyz ?
!
fast3d = ( e1(2) == 0.d0 .and. e1(3) == 0.d0) .and. &
( e2(1) == 0.d0 .and. e2(3) == 0.d0) .and. &
( e3(1) == 0.d0 .and. e3(2) == 0.d0)
!
! are crystal axis aligned along xyz ?
!
fast3d = fast3d .and. &
( at(2,1) == 0.d0 .and. at(3,1) == 0.d0) .and. &
( at(1,2) == 0.d0 .and. at(3,2) == 0.d0) .and. &
( at(1,3) == 0.d0 .and. at(2,3) == 0.d0)
fast3d = fast3d .and. (trim(interpolation) == 'fourier')
!
! Initialise FFT for rho(r) => rho(G) conversion if needed
!
IF (.not. ( iflag == 3 .and. ( output_format == 5 .or. &
output_format == 6 .or. &
fast3d ) ) ) THEN
IF (plot_num==-1) THEN
!
gamma_only=.false.
! nproc_pool=1
!
CALL data_structure ( gamma_only )
CALL allocate_fft()
!
! and rebuild G-vectors in reciprocal space
!
CALL ggen ( gamma_only, at, bg )
!
! here we compute the fourier components of the quantity to plot
!
ELSE
!
IF (gamma_only .and. (trim(interpolation) == 'fourier')) THEN
WRITE(stdout,'(/"BEWARE: plot requiring G-space interpolation",&
&" not implemented for Gamma only!",/, &
&"SOLUTION: restart this calculation with", &
&" emtpy namelist &inputpp")')
CALL errore ('chdens','Not implemented, please read above',1)
ENDIF
!
ENDIF
#if defined(__MPI)
ALLOCATE(aux(dfftp%nnr))
CALL scatter_grid(dfftp, rhor, aux)
psic(:) = cmplx(aux(:), 0.d0,kind=DP)
DEALLOCATE(aux)
ALLOCATE(aux(dfftp%nnr))
CALL scatter_grid(dfftp, rhor, aux)
psic(:) = cmplx(aux(:), 0.d0,kind=DP)
DEALLOCATE(aux)
#else
psic(:) = cmplx(rhor(:), 0.d0,kind=DP)
psic(:) = cmplx(rhor(:), 0.d0,kind=DP)
#endif
CALL fwfft ('Dense', psic, dfftp)
!
! we store the fourier components in the array rhog
!
ALLOCATE (rhog( ngm))
rhog (:) = psic (nl (:) )
!
ENDIF
!
! And now the plot (rhog in G-space, rhor in real space)
!
IF (iflag <= 1) THEN
CALL fwfft ('Dense', psic, dfftp)
!
! we store the fourier components in the array rhog
!
ALLOCATE (rhog( ngm))
rhog (:) = psic (nl (:) )
!
ENDIF
!
! And now the plot (rhog in G-space, rhor in real space)
!
IF (iflag <= 1) THEN
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_1d (nx, m1, x0, e1, ngm, g, rhog, alat, iflag, ounit)
ELSE
CALL plot_1d_bspline (nx, m1, x0, e1, rhor, alat, iflag, ounit)
ENDIF
ELSEIF (iflag == 2) THEN
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_2d (nx, ny, m1, m2, x0, e1, e2, ngm, g, rhog, alat, &
at, nat, tau, atm, ityp, output_format, ounit)
ELSE
CALL plot_2d_bspline (nx, ny, m1, m2, x0, e1, e2, rhor, alat, &
at, nat, tau, atm, ityp, output_format, ounit)
ENDIF
IF (output_format == 2.and.ionode) THEN
WRITE (ounit, '(i4)') nat
WRITE (ounit, '(3f8.4,i3)') ( (tau(ipol,na), ipol=1,3), 1, na=1,nat)
WRITE (ounit, '(f10.6)') celldm (1)
WRITE (ounit, '(3(3f12.6/))') at
ENDIF
ELSEIF (iflag == 3) THEN
IF (output_format == 4.and.ionode) THEN
! gopenmol wants the coordinates in a separate file
IF (fileout /= ' ') THEN
OPEN (unit = ounit+1, file = trim(fileout)//'.xyz', &
form = 'formatted', status = 'unknown')
WRITE( stdout, '(5x,"Writing coordinates to file ",a)') &
trim(fileout)//'.xyz'
ELSE
OPEN (unit = ounit+1, file = 'coord.xyz', &
form = 'formatted', status = 'unknown')
WRITE( stdout, '("Writing coordinates to file coord.xyz")')
ENDIF
ENDIF
IF (output_format == 5.and.ionode) THEN
!
! XCRYSDEN FORMAT
!
CALL xsf_struct (alat, at, nat, tau, atm, ityp, ounit)
CALL xsf_fast_datagrid_3d &
(rhor, dfftp%nr1, dfftp%nr2, dfftp%nr3, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, at, alat, ounit)
ELSEIF (output_format == 6.and.ionode ) THEN
!
! GAUSSIAN CUBE FORMAT
!
IF (TRIM(interpolation) == 'fourier') THEN
CALL write_cubefile (alat, at, bg, nat, tau, atm, ityp, rhor, &
dfftp%nr1, dfftp%nr2, dfftp%nr3, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, ounit)
ELSE
CALL plot_3d_bspline(celldm(1), at, nat, tau, atm, ityp, rhor,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
END IF
ELSEIF (ionode) THEN
!
! GOPENMOL OR XCRYSDEN FORMAT
!
IF (fast3d) THEN
CALL plot_fast (celldm (1), at, nat, tau, atm, ityp, &
dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, dfftp%nr1, dfftp%nr2, dfftp%nr3, rhor, &
bg, m1, m2, m3, x0, e1, e2, e3, output_format, ounit, &
rhotot)
ELSE
IF (nx<=0 .or. ny <=0 .or. nz <=0) &
CALL errore("chdens","nx,ny,nz, required",1)
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_3d (celldm (1), at, nat, tau, atm, ityp, ngm, g, rhog,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
ELSE
CALL plot_3d_bspline(celldm(1), at, nat, tau, atm, ityp, rhor,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
ENDIF
!
ENDIF
ENDIF
ELSEIF (iflag == 4) THEN
radius = radius / alat
CALL plot_2ds (nx, ny, radius, ngm, g, rhog, output_format, ounit)
ELSE
CALL errore ('chdens', 'wrong iflag', 1)
ENDIF
!
WRITE(stdout, '(5x,"Plot Type: ",a," Output format: ",a)') &
plotname(iflag), formatname(output_format)
!
IF (allocated(rhog)) DEALLOCATE(rhog)
DEALLOCATE(rhor)
ENDDO
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_1d (nx, m1, x0, e1, ngm, g, rhog, alat, iflag, ounit)
ELSE
CALL plot_1d_bspline (nx, m1, x0, e1, rhor, alat, iflag, ounit)
ENDIF
ELSEIF (iflag == 2) THEN
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_2d (nx, ny, m1, m2, x0, e1, e2, ngm, g, rhog, alat, &
at, nat, tau, atm, ityp, output_format, ounit)
ELSE
CALL plot_2d_bspline (nx, ny, m1, m2, x0, e1, e2, rhor, alat, &
at, nat, tau, atm, ityp, output_format, ounit)
ENDIF
IF (output_format == 2.and.ionode) THEN
WRITE (ounit, '(i4)') nat
WRITE (ounit, '(3f8.4,i3)') ( (tau(ipol,na), ipol=1,3), 1, na=1,nat)
WRITE (ounit, '(f10.6)') celldm (1)
WRITE (ounit, '(3(3f12.6/))') at
ENDIF
ELSEIF (iflag == 3) THEN
IF (output_format == 4.and.ionode) THEN
! gopenmol wants the coordinates in a separate file
IF (fileout /= ' ') THEN
OPEN (unit = ounit+1, file = trim(fileout)//'.xyz', &
form = 'formatted', status = 'unknown')
WRITE( stdout, '(5x,"Writing coordinates to file ",a)') &
trim(fileout)//'.xyz'
ELSE
OPEN (unit = ounit+1, file = 'coord.xyz', &
form = 'formatted', status = 'unknown')
WRITE( stdout, '("Writing coordinates to file coord.xyz")')
ENDIF
ENDIF
IF (output_format == 5.and.ionode) THEN
!
! XCRYSDEN FORMAT
!
CALL xsf_struct (alat, at, nat, tau, atm, ityp, ounit)
CALL xsf_fast_datagrid_3d &
(rhor, dfftp%nr1, dfftp%nr2, dfftp%nr3, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, at, alat, ounit)
ELSEIF (output_format == 6.and.ionode ) THEN
!
! GAUSSIAN CUBE FORMAT
!
IF (TRIM(interpolation) == 'fourier') THEN
CALL write_cubefile (alat, at, bg, nat, tau, atm, ityp, rhor, &
dfftp%nr1, dfftp%nr2, dfftp%nr3, dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, ounit)
ELSE
CALL plot_3d_bspline(celldm(1), at, nat, tau, atm, ityp, rhor,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
END IF
ELSEIF (ionode) THEN
!
! GOPENMOL OR XCRYSDEN FORMAT
!
IF (fast3d) THEN
CALL plot_fast (celldm (1), at, nat, tau, atm, ityp, &
dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, dfftp%nr1, dfftp%nr2, dfftp%nr3, rhor, &
bg, m1, m2, m3, x0, e1, e2, e3, output_format, ounit, &
rhotot)
ELSE
IF (nx<=0 .or. ny <=0 .or. nz <=0) &
CALL errore("chdens","nx,ny,nz, required",1)
IF (TRIM(interpolation) == 'fourier') THEN
CALL plot_3d (celldm (1), at, nat, tau, atm, ityp, ngm, g, rhog,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
ELSE
CALL plot_3d_bspline(celldm(1), at, nat, tau, atm, ityp, rhor,&
nx, ny, nz, m1, m2, m3, x0, e1, e2, e3, output_format, &
ounit, rhotot)
ENDIF
!
ENDIF
ENDIF
ELSEIF (iflag == 4) THEN
radius = radius / alat
CALL plot_2ds (nx, ny, radius, ngm, g, rhog, output_format, ounit)
ELSE
CALL errore ('chdens', 'wrong iflag', 1)
ENDIF
!
WRITE(stdout, '(5x,"Plot Type: ",a," Output format: ",a)') &
plotname(iflag), formatname(output_format)
!
IF (allocated(rhog)) DEALLOCATE(rhog)
DEALLOCATE(rhor)
DEALLOCATE(tau)
DEALLOCATE(ityp)
@ -1449,3 +1471,5 @@ SUBROUTINE isostm_plot(rhor, nr1x, nr2x, nr3x, &
DEALLOCATE(image)
END SUBROUTINE isostm_plot
END MODULE chdens_module

41
PP/src/local_dos.f90
View File

@ -35,7 +35,7 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
nkstot, ngk, igk_k
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE scf, ONLY : rho
USE symme, ONLY : sym_rho, sym_rho_init
USE symme, ONLY : sym_rho, sym_rho_init, sym_rho_deallocate
USE uspp, ONLY : nkb, vkb, becsum, nhtol, nhtoj, indv
USE uspp_param, ONLY : upf, nh, nhm
USE wavefunctions_module, ONLY : evc, psic, psic_nc
@ -56,23 +56,24 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
!
INTEGER, INTENT(in) :: iflag, kpoint, kband, spin_component
LOGICAL, INTENT(in) :: lsign
real(DP), INTENT(in) :: emin, emax
REAL(DP), INTENT(in) :: emin, emax
!
real(DP), INTENT(out) :: dos (dfftp%nnr)
REAL(DP), INTENT(out) :: dos (dfftp%nnr)
!
! local variables
!
INTEGER :: npw, ikb, jkb, ijkb0, ih, jh, kh, na, ijh, np
! counters for US PPs
INTEGER :: ir, is, ig, ibnd, ik, irm, isup, isdw, ipol, kkb, is1, is2
INTEGER :: npw, ikb, jkb, ijkb0, ih, jh, kh, na, ijh, np
! counters
real(DP) :: w, w1, modulus
real(DP), ALLOCATABLE :: rbecp(:,:), segno(:), maxmod(:)
INTEGER :: ir, is, ig, ibnd, ik, irm, isup, isdw, ipol, kkb, is1, is2
REAL(DP) :: w, w1, modulus, wg_max
REAL(DP), ALLOCATABLE :: rbecp(:,:), segno(:), maxmod(:)
COMPLEX(DP), ALLOCATABLE :: becp(:,:), &
becp_nc(:,:,:), be1(:,:), be2(:,:)
INTEGER :: who_calculate, iproc
COMPLEX(DP) :: phase
real(DP), EXTERNAL :: w0gauss, w1gauss
REAL(DP), EXTERNAL :: w0gauss, w1gauss
LOGICAL :: i_am_the_pool
INTEGER :: which_pool, kpoint_pool
!
@ -112,16 +113,16 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
! calculate the correct weights
!
IF (iflag /= 0.and. iflag /=3 .and. .not.lgauss) CALL errore ('local_dos', &
'gaussian broadening needed', 1)
'gaussian broadening needed', 1)
IF (iflag == 2 .and. ngauss /= -99) CALL errore ('local_dos', &
' beware: not using Fermi-Dirac function ', - ngauss)
' beware: not using Fermi-Dirac function ', - ngauss)
DO ik = 1, nks
DO ibnd = 1, nbnd
IF (iflag == 0) THEN
wg (ibnd, ik) = 0.d0
ELSEIF (iflag == 1) THEN
wg (ibnd, ik) = wk (ik) * w0gauss ( (ef - et (ibnd, ik) ) &
/ degauss, ngauss) / degauss
! Local density of states at energy emin with broadening emax
wg(ibnd,ik) = wk(ik) * w0gauss((emin - et(ibnd, ik))/emax, ngauss) / emax
ELSEIF (iflag == 2) THEN
wg (ibnd, ik) = - wk (ik) * w1gauss ( (ef - et (ibnd, ik) ) &
/ degauss, ngauss)
@ -136,6 +137,7 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
ENDIF
ENDDO
ENDDO
wg_max = MAXVAL(wg(:,:))
IF ( iflag == 0 .and. npool > 1 ) THEN
CALL xk_pool( kpoint, nkstot, kpoint_pool, which_pool )
@ -170,7 +172,9 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
! here we compute the density of states
!
DO ibnd = 1, nbnd
IF (ibnd == kband .or. iflag /= 0) THEN
! Neglect summands with relative weights below machine epsilon
IF ( wg(ibnd, ik) > epsilon(0.0_DP) * wg_max .and. &
(ibnd == kband .or. iflag /= 0)) THEN
IF (noncolin) THEN
psic_nc = (0.d0,0.d0)
DO ig = 1, npw
@ -352,9 +356,10 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
ENDIF
ENDDO
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO ! loop over bands
ENDIF
ENDDO ! loop over k-points
IF (gamma_only) THEN
DEALLOCATE(rbecp)
ELSE
@ -412,7 +417,7 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
!
! symmetrization of the local dos
!
CALL sym_rho_init ( gamma_only )
CALL sym_rho_init (gamma_only )
!
psic(:) = cmplx ( dos(:), 0.0_dp, kind=dp)
CALL fwfft ('Dense', psic, dfftp)
@ -425,6 +430,8 @@ SUBROUTINE local_dos (iflag, lsign, kpoint, kband, spin_component, &
CALL invfft ('Dense', psic, dfftp)
dos(:) = dble(psic(:))
!
CALL sym_rho_deallocate()
!
RETURN
END SUBROUTINE local_dos

48
PP/src/make.depend
View File

@ -110,33 +110,34 @@ bgw2pw.o : ../../PW/src/scf_mod.o
bgw2pw.o : ../../PW/src/symm_base.o
bgw2pw.o : ../../iotk/src/iotk_module.o
cft.o : ../../Modules/kind.o
chdens.o : ../../FFTXlib/fft_interfaces.o
chdens.o : ../../FFTXlib/fft_types.o
chdens.o : ../../FFTXlib/scatter_mod.o
chdens.o : ../../Modules/cell_base.o
chdens.o : ../../Modules/constants.o
chdens.o : ../../Modules/control_flags.o
chdens.o : ../../Modules/fft_base.o
chdens.o : ../../Modules/gvecw.o
chdens.o : ../../Modules/io_files.o
chdens.o : ../../Modules/io_global.o
chdens.o : ../../Modules/ions_base.o
chdens.o : ../../Modules/kind.o
chdens.o : ../../Modules/mp.o
chdens.o : ../../Modules/mp_bands.o
chdens.o : ../../Modules/mp_pools.o
chdens.o : ../../Modules/mp_world.o
chdens.o : ../../Modules/parameters.o
chdens.o : ../../Modules/recvec.o
chdens.o : ../../Modules/recvec_subs.o
chdens.o : ../../Modules/run_info.o
chdens.o : ../../Modules/wavefunctions.o
chdens.o : ../../PW/src/pwcom.o
chdens_bspline.o : ../../Modules/bspline.o
chdens_bspline.o : ../../Modules/cell_base.o
chdens_bspline.o : ../../Modules/fft_base.o
chdens_bspline.o : ../../Modules/io_global.o
chdens_bspline.o : ../../Modules/kind.o
chdens_bspline.o : chdens_module.o
chdens_module.o : ../../FFTXlib/fft_interfaces.o
chdens_module.o : ../../FFTXlib/fft_types.o
chdens_module.o : ../../FFTXlib/scatter_mod.o
chdens_module.o : ../../Modules/cell_base.o
chdens_module.o : ../../Modules/constants.o
chdens_module.o : ../../Modules/control_flags.o
chdens_module.o : ../../Modules/fft_base.o
chdens_module.o : ../../Modules/gvecw.o
chdens_module.o : ../../Modules/io_files.o
chdens_module.o : ../../Modules/io_global.o
chdens_module.o : ../../Modules/ions_base.o
chdens_module.o : ../../Modules/kind.o
chdens_module.o : ../../Modules/mp.o
chdens_module.o : ../../Modules/mp_bands.o
chdens_module.o : ../../Modules/mp_pools.o
chdens_module.o : ../../Modules/mp_world.o
chdens_module.o : ../../Modules/parameters.o
chdens_module.o : ../../Modules/recvec.o
chdens_module.o : ../../Modules/recvec_subs.o
chdens_module.o : ../../Modules/run_info.o
chdens_module.o : ../../Modules/wavefunctions.o
chdens_module.o : ../../PW/src/pwcom.o
compute_ppsi.o : ../../Modules/becmod.o
compute_ppsi.o : ../../Modules/cell_base.o
compute_ppsi.o : ../../Modules/io_global.o
@ -163,6 +164,7 @@ compute_sigma_avg.o : ../../PW/src/scf_mod.o
cube.o : ../../Modules/cell_base.o
cube.o : ../../Modules/io_global.o
cube.o : ../../Modules/kind.o
cube.o : ../../Modules/run_info.o
d_matrix_nc.o : ../../Modules/invmat.o
d_matrix_nc.o : ../../Modules/kind.o
d_matrix_nc.o : ../../Modules/random_numbers.o
@ -429,9 +431,11 @@ postproc.o : ../../Modules/mp.o
postproc.o : ../../Modules/mp_global.o
postproc.o : ../../Modules/mp_world.o
postproc.o : ../../Modules/noncol.o
postproc.o : ../../Modules/parameters.o
postproc.o : ../../Modules/paw_variables.o
postproc.o : ../../Modules/recvec.o
postproc.o : ../../PW/src/pwcom.o
postproc.o : chdens_module.o
projections_mod.o : ../../Modules/ions_base.o
projections_mod.o : ../../Modules/kind.o
projections_mod.o : ../../Modules/noncol.o

212
PP/src/postproc.f90
View File

@ -7,50 +7,10 @@
!
!
!-----------------------------------------------------------------------
PROGRAM pp
!-----------------------------------------------------------------------
!
! Program for data analysis and plotting. The two basic steps are:
! 1) read the output file produced by pw.x, extract and calculate
! the desired quantity (rho, V, ...)
! 2) write the desired quantity to file in a suitable format for
! various types of plotting and various plotting programs
! The two steps can be performed independently. Intermediate data
! can be saved to file in step 1 and read from file in step 2.
!
! DESCRIPTION of the INPUT : see file Doc/INPUT_PP.*
!
USE io_global, ONLY : ionode
USE mp_global, ONLY : mp_startup
USE environment,ONLY : environment_start, environment_end
!
IMPLICIT NONE
!
CHARACTER(len=256) :: filplot
INTEGER :: plot_num
!
! initialise environment
!
#if defined(__MPI)
CALL mp_startup ( )
#endif
CALL environment_start ( 'POST-PROC' )
!
IF ( ionode ) CALL input_from_file ( )
!
CALL extract (filplot, plot_num)
!
CALL chdens (filplot, plot_num)
!
CALL environment_end ( 'POST-PROC' )
!
CALL stop_pp()
!
END PROGRAM pp
!
MODULE pp_module
CONTAINS
!-----------------------------------------------------------------------
SUBROUTINE extract (filplot,plot_num)
SUBROUTINE extract (plot_files,plot_num)
!-----------------------------------------------------------------------
!
! This subroutine reads the data for the output file produced by pw.x
@ -65,7 +25,7 @@ SUBROUTINE extract (filplot,plot_num)
USE ions_base, ONLY : nat, ntyp=>nsp, ityp, tau
USE gvect
USE fft_base, ONLY : dfftp
USE klist, ONLY : two_fermi_energies
USE klist, ONLY : two_fermi_energies, degauss
USE vlocal, ONLY : strf
USE io_files, ONLY : tmp_dir, prefix
USE io_global, ONLY : ionode, ionode_id
@ -76,25 +36,35 @@ SUBROUTINE extract (filplot,plot_num)
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE constants, ONLY : rytoev
USE parameters, ONLY : npk
USE io_global, ONLY : stdout
IMPLICIT NONE
!
CHARACTER(LEN=256), EXTERNAL :: trimcheck
!
CHARACTER(len=256), INTENT(out) :: filplot
CHARACTER(len=256), DIMENSION(:), ALLOCATABLE, INTENT(out) :: plot_files
INTEGER, INTENT(out) :: plot_num
INTEGER :: kpoint, kband, spin_component, ios
CHARACTER (len=2), DIMENSION(0:3) :: spin_desc = &
(/ ' ', '_X', '_Y', '_Z' /)
INTEGER :: kpoint(2), kband(2), spin_component(3), ios
LOGICAL :: lsign, needwf
REAL(DP) :: emin, emax, sample_bias, z, dz, epsilon
REAL(DP) :: degauss_ldos, delta_e
CHARACTER(len=256) :: filplot
INTEGER :: plot_nkpt, plot_nbnd, plot_nspin, nplots
INTEGER :: iplot, ikpt, ibnd, ispin
! directory for temporary files
CHARACTER(len=256) :: outdir
NAMELIST / inputpp / outdir, prefix, plot_num, sample_bias, &
spin_component, z, dz, emin, emax, kpoint, kband, &
filplot, lsign, epsilon
spin_component, z, dz, emin, emax, delta_e, degauss_ldos, kpoint, kband, &
filplot, lsign, epsilon
!
! set default values for variables in namelist
!
@ -103,6 +73,8 @@ SUBROUTINE extract (filplot,plot_num)
IF ( trim( outdir ) == ' ' ) outdir = './'
filplot = 'tmp.pp'
plot_num = -1
kpoint(2) = 0
kband(2) = 0
spin_component = 0
sample_bias = 0.01d0
z = 1.d0
@ -110,7 +82,9 @@ SUBROUTINE extract (filplot,plot_num)
lsign=.false.
emin = -999.0d0
emax = +999.0d0
delta_e=0.1d0
epsilon=1.d0
degauss_ldos=-999.0d0
!
ios = 0
!
@ -139,6 +113,8 @@ SUBROUTINE extract (filplot,plot_num)
CALL mp_bcast( dz, ionode_id, world_comm )
CALL mp_bcast( emin, ionode_id, world_comm )
CALL mp_bcast( emax, ionode_id, world_comm )
CALL mp_bcast( degauss_ldos, ionode_id, world_comm )
CALL mp_bcast( delta_e, ionode_id, world_comm )
CALL mp_bcast( kband, ionode_id, world_comm )
CALL mp_bcast( kpoint, ionode_id, world_comm )
CALL mp_bcast( filplot, ionode_id, world_comm )
@ -153,13 +129,13 @@ SUBROUTINE extract (filplot,plot_num)
'Wrong plot_num', abs (plot_num) )
IF (plot_num == 7 .or. plot_num == 13 .or. plot_num==18) THEN
IF (spin_component < 0 .or. spin_component > 3) CALL errore &
IF (spin_component(1) < 0 .or. spin_component(1) > 3) CALL errore &
('postproc', 'wrong spin_component', 1)
ELSEIF (plot_num == 10) THEN
IF (spin_component < 0 .or. spin_component > 2) CALL errore &
IF (spin_component(1) < 0 .or. spin_component(1) > 2) CALL errore &
('postproc', 'wrong spin_component', 2)
ELSE
IF (spin_component < 0 ) CALL errore &
IF (spin_component(1) < 0 ) CALL errore &
('postproc', 'wrong spin_component', 3)
ENDIF
!
@ -182,18 +158,128 @@ SUBROUTINE extract (filplot,plot_num)
CALL errore('postproc',&
'Post-processing with constrained magnetization is not available yet',1)
!
! The following line sets emax to its default value if not set
! It is done here because Ef must be read from file
!
IF (emax == +999.0d0) emax = ef
! Set default values for emin, emax, degauss_ldos
! Done here because ef, degauss must be read from file
IF (emin > emax) CALL errore('postproc','emin > emax',0)
IF (plot_num == 10) THEN
emin = emin / rytoev
emax = emax / rytoev
IF (emax == +999.0d0) emax = ef * rytoev
ELSEIF (plot_num == 3) THEN
IF (emin == -999.0d0) emin = ef * rytoev
IF (emax == +999.0d0) emax = ef * rytoev
IF (degauss_ldos == -999.0d0) THEN
WRITE(stdout, &
'(/5x,"degauss_ldos not set, defaults to degauss = ",f6.4, " eV")') &
degauss * rytoev
degauss_ldos = degauss * rytoev
ENDIF
ENDIF
! transforming all back to Ry units
emin = emin / rytoev
emax = emax / rytoev
delta_e = delta_e / rytoev
degauss_ldos = degauss_ldos / rytoev
! Number of output files depends on input
nplots = 1
IF (plot_num == 3) THEN
nplots=(emax-emin)/delta_e + 1
ELSEIF (plot_num == 7) THEN
IF (kpoint(2) == 0) kpoint(2) = kpoint(1)
plot_nkpt = kpoint(2) - kpoint(1) + 1
IF (kband(2) == 0) kband(2) = kband(1)
plot_nbnd = kband(2) - kband(1) + 1
IF (spin_component(2) == 0) spin_component(2) = spin_component(1)
plot_nspin = spin_component(2) - spin_component(1) + 1
nplots = plot_nbnd * plot_nkpt * plot_nspin
ENDIF
ALLOCATE( plot_files(nplots) )
plot_files(1) = filplot
!
! First handle plot_nums with multiple calls to punch_plot
!
IF (nplots > 1 .AND. plot_num == 3) THEN
! Local density of states on energy grid of spacing delta_e within [emin, emax]
DO iplot=1,nplots
WRITE(plot_files(iplot),'(A, I0.3)') TRIM(filplot), iplot
CALL punch_plot (TRIM(plot_files(iplot)), plot_num, sample_bias, z, dz, &
emin, degauss_ldos, kpoint, kband, spin_component, lsign, epsilon)
emin=emin+delta_e
ENDDO
ELSEIF (nplots > 1 .AND. plot_num == 7) THEN
! Plot multiple KS orbitals in one go
iplot = 1
DO ikpt=kpoint(1), kpoint(2)
DO ibnd=kband(1), kband(2)
DO ispin=spin_component(1), spin_component(2)
WRITE(plot_files(iplot),"(A,A,I0.3,A,I0.3,A)") &
TRIM(filplot), "_K", ikpt, "_B", ibnd, TRIM(spin_desc(ispin))
CALL punch_plot (TRIM(plot_files(iplot)), plot_num, sample_bias, z, dz, &
emin, emax, ikpt, ibnd, ispin, lsign, epsilon)
iplot = iplot + 1
ENDDO
ENDDO
ENDDO
ELSE
! Single call to punch_plot
IF (plot_num == 3) THEN
CALL punch_plot (filplot, plot_num, sample_bias, z, dz, &
emin, degauss_ldos, kpoint, kband, spin_component, lsign, epsilon)
ELSE
CALL punch_plot (filplot, plot_num, sample_bias, z, dz, &
emin, emax, kpoint, kband, spin_component, lsign, epsilon)
ENDIF
ENDIF
!
! Now do whatever you want
!
CALL punch_plot (filplot, plot_num, sample_bias, z, dz, &
emin, emax, kpoint, kband, spin_component, lsign, epsilon)
!
END SUBROUTINE extract
END MODULE pp_module
!
!-----------------------------------------------------------------------
PROGRAM pp
!-----------------------------------------------------------------------
!
! Program for data analysis and plotting. The two basic steps are:
! 1) read the output file produced by pw.x, extract and calculate
! the desired quantity (rho, V, ...)
! 2) write the desired quantity to file in a suitable format for
! various types of plotting and various plotting programs
! The two steps can be performed independently. Intermediate data
! can be saved to file in step 1 and read from file in step 2.
!
! DESCRIPTION of the INPUT : see file Doc/INPUT_PP.*
!
USE io_global, ONLY : ionode
USE mp_global, ONLY : mp_startup
USE environment,ONLY : environment_start, environment_end
USE chdens_module, ONLY : chdens
USE pp_module, ONLY : extract
!
IMPLICIT NONE
!
!CHARACTER(len=256) :: filplot
CHARACTER(len=256), DIMENSION(:), ALLOCATABLE :: plot_files
INTEGER :: plot_num
!
! initialise environment
!
#if defined(__MPI)
CALL mp_startup ( )
#endif
CALL environment_start ( 'POST-PROC' )
!
IF ( ionode ) CALL input_from_file ( )
!
CALL extract (plot_files, plot_num)
!
CALL chdens (plot_files, plot_num)
!
CALL environment_end ( 'POST-PROC' )
!
CALL stop_pp()
!
END PROGRAM pp

27
PP/src/punch_plot.f90
View File

@ -40,11 +40,12 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
USE paw_postproc, ONLY : PAW_make_ae_charge
IMPLICIT NONE
CHARACTER(len=*) :: filplot
INTEGER :: kpoint, kband, spin_component, plot_num
LOGICAL :: lsign
REAL(DP) :: sample_bias, dummy
REAL(DP) :: emin, emax, z, dz, charge, epsilon
CHARACTER(len=*), INTENT(IN) :: filplot
INTEGER, INTENT(IN) :: plot_num, kpoint, kband, spin_component
LOGICAL, INTENT(IN) :: lsign
REAL(DP), INTENT(IN) :: sample_bias, z, dz, &
emin, emax, epsilon
REAL(DP) :: dummy, charge
INTEGER :: is, ipol, istates
#if defined(__MPI)
! auxiliary vector (parallel case)
@ -60,10 +61,13 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
#endif
WRITE( stdout, '(/5x,"Calling punch_plot, plot_num = ",i3)') plot_num
IF (plot_num == 7 ) &
WRITE( stdout, '(/5x,"Plotting k_point = ",i3," band =", i3 )') &
IF (plot_num == 3 ) &
WRITE(stdout, '(/5x,"Energy =", f10.5, " eV, broadening =", f10.5, "eV" )') &
emin * rytoev, emax * rytoev
IF (plot_num == 7) &
WRITE( stdout, '(/5x,"Plotting k_point = ",i3," band =", i3 )') &
kpoint, kband
IF (plot_num == 7 .and. noncolin .and. spin_component /= 0 ) &
IF ((plot_num == 7) .and. noncolin .and. spin_component /= 0 ) &
WRITE( stdout, '(/5x,"Plotting spin magnetization ipol = ",i3)') &
spin_component
!
@ -118,10 +122,12 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
ELSEIF (plot_num == 3) THEN
!
! The local density of states at e_fermi on output
! The local density of states at emin, with broadening emax
!
WRITE (title, '(" Energy = ",f8.4," eV, ", "broadening = ",f8.4," eV")') &
emin * rytoev, emax * rytoev
IF (noncolin) CALL errore('punch_plot','not implemented yet',1)
CALL local_dos (1, lsign, kpoint, kband, spin_component, emin, emax, raux)
CALL local_dos(1, lsign, kpoint, kband, spin_component, emin, emax, raux)
ELSEIF (plot_num == 4) THEN
!
@ -153,6 +159,7 @@ SUBROUTINE punch_plot (filplot, plot_num, sample_bias, z, dz, &
ENDIF
ELSEIF (plot_num == 7) THEN
WRITE (title, '("k_point ",i4,", band ",i4)') kpoint ,kband
IF (noncolin) THEN
IF (spin_component==0) THEN