1!
2! Copyright (C) 2001-2009 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!-----------------------------------------------------------------------
10MODULE pp_module
11CONTAINS
12!-----------------------------------------------------------------------
13SUBROUTINE extract (plot_files,plot_num)
14  !-----------------------------------------------------------------------
15  !
16  !    Reads data produced by pw.x, computes the desired quantity (rho, V, ...)
17  !    and writes it to a file (or multiple files) for further processing or
18  !    plotting
19  !
20  !    On return, plot_files contains a list of all written files.
21  !
22  !    DESCRIPTION of the INPUT: see file Doc/INPUT_PP
23  !
24  USE kinds,     ONLY : DP
25  USE cell_base, ONLY : bg
26  USE ener,      ONLY : ef
27  USE ions_base, ONLY : nat, ntyp=>nsp, ityp, tau
28  USE gvect
29  USE fft_base,  ONLY : dfftp
30  USE klist,     ONLY : two_fermi_energies, degauss
31  USE vlocal,    ONLY : strf
32  USE io_files,  ONLY : tmp_dir, prefix
33  USE io_global, ONLY : ionode, ionode_id
34  USE noncollin_module, ONLY : i_cons
35  USE paw_variables, ONLY : okpaw
36  USE mp,        ONLY : mp_bcast
37  USE mp_images, ONLY : intra_image_comm
38  USE constants, ONLY : rytoev
39  USE parameters, ONLY : npk
40  USE io_global, ONLY : stdout
41  !
42  IMPLICIT NONE
43  !
44  CHARACTER(LEN=256), EXTERNAL :: trimcheck
45  !
46  CHARACTER(len=256), DIMENSION(:), ALLOCATABLE, INTENT(out) :: plot_files
47  INTEGER, INTENT(out) :: plot_num
48
49  CHARACTER (len=2), DIMENSION(0:3) :: spin_desc = &
50       (/ '  ', '_X', '_Y', '_Z' /)
51
52  INTEGER :: kpoint(2), kband(2), spin_component(3), ios
53  LOGICAL :: lsign, needwf, dummy
54
55  REAL(DP) :: emin, emax, sample_bias, z, dz
56
57  REAL(DP) :: degauss_ldos, delta_e
58  CHARACTER(len=256) :: filplot
59  INTEGER :: plot_nkpt, plot_nbnd, plot_nspin, nplots
60  INTEGER :: iplot, ikpt, ibnd, ispin
61
62  ! directory for temporary files
63  CHARACTER(len=256) :: outdir
64
65  NAMELIST / inputpp / outdir, prefix, plot_num, sample_bias, &
66      spin_component, z, dz, emin, emax, delta_e, degauss_ldos, kpoint, kband, &
67      filplot, lsign
68  !
69  !   set default values for variables in namelist
70  !
71  prefix = 'pwscf'
72  CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir )
73  IF ( trim( outdir ) == ' ' ) outdir = './'
74  filplot = 'tmp.pp'
75  plot_num = -1
76  kpoint(2) = 0
77  kband(2) = 0
78  spin_component = 0
79  sample_bias = 0.01d0
80  z = 1.d0
81  dz = 0.05d0
82  lsign=.false.
83  emin = -999.0d0
84  emax = +999.0d0
85  delta_e=0.1d0
86  degauss_ldos=-999.0d0
87  !
88  ios = 0
89  !
90  IF ( ionode )  THEN
91     !
92     !     reading the namelist inputpp
93     !
94     READ (5, inputpp, iostat = ios)
95     !
96     tmp_dir = trimcheck ( outdir )
97     !
98  ENDIF
99  !
100  CALL mp_bcast (ios, ionode_id, intra_image_comm)
101  !
102  IF ( ios /= 0) CALL errore ('postproc', 'reading inputpp namelist', abs(ios))
103  !
104  ! ... Broadcast variables
105  !
106  CALL mp_bcast( tmp_dir, ionode_id, intra_image_comm )
107  CALL mp_bcast( prefix, ionode_id, intra_image_comm )
108  CALL mp_bcast( plot_num, ionode_id, intra_image_comm )
109  CALL mp_bcast( sample_bias, ionode_id, intra_image_comm )
110  CALL mp_bcast( spin_component, ionode_id, intra_image_comm )
111  CALL mp_bcast( z, ionode_id, intra_image_comm )
112  CALL mp_bcast( dz, ionode_id, intra_image_comm )
113  CALL mp_bcast( emin, ionode_id, intra_image_comm )
114  CALL mp_bcast( emax, ionode_id, intra_image_comm )
115  CALL mp_bcast( degauss_ldos, ionode_id, intra_image_comm )
116  CALL mp_bcast( delta_e, ionode_id, intra_image_comm )
117  CALL mp_bcast( kband, ionode_id, intra_image_comm )
118  CALL mp_bcast( kpoint, ionode_id, intra_image_comm )
119  CALL mp_bcast( filplot, ionode_id, intra_image_comm )
120  CALL mp_bcast( lsign, ionode_id, intra_image_comm )
121  !
122  ! no task specified: do nothing and return
123  !
124  IF (plot_num == -1) THEN
125     ALLOCATE( plot_files(0) )
126     RETURN
127  ENDIF
128  !
129  IF (plot_num < 0 .or. plot_num > 22) CALL errore ('postproc', &
130          'Wrong plot_num', abs (plot_num) )
131
132  IF (plot_num == 7 .or. plot_num == 13 .or. plot_num==18) THEN
133     IF  (spin_component(1) < 0 .or. spin_component(1) > 3) CALL errore &
134          ('postproc', 'wrong spin_component', 1)
135  ELSEIF (plot_num == 10) THEN
136     IF  (spin_component(1) < 0 .or. spin_component(1) > 2) CALL errore &
137          ('postproc', 'wrong spin_component', 2)
138  ELSE
139     IF (spin_component(1) < 0 ) CALL errore &
140         ('postproc', 'wrong spin_component', 3)
141  ENDIF
142  !
143  !   Read xml file, allocate and initialize general variables
144  !   If needed, allocate and initialize wavefunction-related variables
145  !
146  needwf=(plot_num==3).or.(plot_num==4).or.(plot_num==5).or.(plot_num==7).or. &
147         (plot_num==8).or.(plot_num==10)
148  CALL read_file_new ( needwf )
149  !
150  IF ( ( two_fermi_energies .or. i_cons /= 0) .and. &
151       ( plot_num==3 .or. plot_num==4 .or. plot_num==5 ) ) &
152     CALL errore('postproc',&
153     'Post-processing with constrained magnetization is not available yet',1)
154  !
155  ! Set default values for emin, emax, degauss_ldos
156  ! Done here because ef, degauss must be read from file
157  IF (emin > emax) CALL errore('postproc','emin > emax',0)
158  IF (plot_num == 10) THEN
159      IF (emax == +999.0d0) emax = ef * rytoev
160  ELSEIF (plot_num == 3) THEN
161      IF (emin == -999.0d0) emin = ef * rytoev
162      IF (emax == +999.0d0) emax = ef * rytoev
163      IF (degauss_ldos == -999.0d0) THEN
164          WRITE(stdout, &
165              '(/5x,"degauss_ldos not set, defaults to degauss = ",f6.4, " eV")') &
166             degauss * rytoev
167          degauss_ldos = degauss * rytoev
168      ENDIF
169  ENDIF
170  ! transforming all back to Ry units
171  emin = emin / rytoev
172  emax = emax / rytoev
173  delta_e = delta_e / rytoev
174  degauss_ldos = degauss_ldos / rytoev
175
176  ! Number of output files depends on input
177  nplots = 1
178  IF (plot_num == 3) THEN
179     nplots=(emax-emin)/delta_e + 1
180  ELSEIF (plot_num == 7) THEN
181      IF (kpoint(2) == 0)  kpoint(2) = kpoint(1)
182      plot_nkpt = kpoint(2) - kpoint(1) + 1
183      IF (kband(2) == 0)  kband(2) = kband(1)
184      plot_nbnd = kband(2) - kband(1) + 1
185      IF (spin_component(2) == 0)  spin_component(2) = spin_component(1)
186      plot_nspin = spin_component(2) - spin_component(1) + 1
187
188      nplots = plot_nbnd * plot_nkpt * plot_nspin
189  ENDIF
190  ALLOCATE( plot_files(nplots) )
191  plot_files(1) = filplot
192
193  !
194  ! First handle plot_nums with multiple calls to punch_plot
195  !
196  IF (nplots > 1 .AND. plot_num == 3) THEN
197  ! Local density of states on energy grid of spacing delta_e within [emin, emax]
198    DO iplot=1,nplots
199      WRITE(plot_files(iplot),'(A, I0.3)') TRIM(filplot), iplot
200      CALL punch_plot (TRIM(plot_files(iplot)), plot_num, sample_bias, z, dz, &
201        emin, degauss_ldos, kpoint, kband, spin_component, lsign)
202      emin=emin+delta_e
203    ENDDO
204  ELSEIF (nplots > 1 .AND. plot_num == 7) THEN
205  ! Plot multiple KS orbitals in one go
206    iplot = 1
207    DO ikpt=kpoint(1), kpoint(2)
208      DO ibnd=kband(1), kband(2)
209        DO ispin=spin_component(1), spin_component(2)
210          WRITE(plot_files(iplot),"(A,A,I0.3,A,I0.3,A)") &
211            TRIM(filplot), "_K", ikpt, "_B", ibnd, TRIM(spin_desc(ispin))
212          CALL punch_plot (TRIM(plot_files(iplot)), plot_num, sample_bias, z, dz, &
213            emin, emax, ikpt, ibnd, ispin, lsign)
214          iplot = iplot + 1
215        ENDDO
216      ENDDO
217    ENDDO
218
219  ELSE
220  ! Single call to punch_plot
221    IF (plot_num == 3) THEN
222       CALL punch_plot (filplot, plot_num, sample_bias, z, dz, &
223           emin, degauss_ldos, kpoint, kband, spin_component, lsign)
224     ELSE
225       CALL punch_plot (filplot, plot_num, sample_bias, z, dz, &
226          emin, emax, kpoint, kband, spin_component, lsign)
227     ENDIF
228
229  ENDIF
230  !
231END SUBROUTINE extract
232
233END MODULE pp_module
234!
235!-----------------------------------------------------------------------
236PROGRAM pp
237  !-----------------------------------------------------------------------
238  !
239  !    Program for data analysis and plotting. The two basic steps are:
240  !    1) read the output file produced by pw.x, extract and calculate
241  !       the desired quantity (rho, V, ...)
242  !    2) write the desired quantity to file in a suitable format for
243  !       various types of plotting and various plotting programs
244  !    The two steps can be performed independently. Intermediate data
245  !    can be saved to file in step 1 and read from file in step 2.
246  !
247  !    DESCRIPTION of the INPUT : see file Doc/INPUT_PP.*
248  !
249  USE io_global,  ONLY : ionode
250  USE mp_global,  ONLY : mp_startup
251  USE environment,ONLY : environment_start, environment_end
252  USE chdens_module, ONLY : chdens
253  USE pp_module, ONLY : extract
254
255  !
256  IMPLICIT NONE
257  !
258  !CHARACTER(len=256) :: filplot
259  CHARACTER(len=256), DIMENSION(:), ALLOCATABLE :: plot_files
260  INTEGER :: plot_num
261  !
262  ! initialise environment
263  !
264#if defined(__MPI)
265  CALL mp_startup ( )
266#endif
267  CALL environment_start ( 'POST-PROC' )
268  !
269  IF ( ionode )  CALL input_from_file ( )
270  !
271  CALL extract (plot_files, plot_num)
272  !
273  CALL chdens (plot_files, plot_num)
274  !
275  CALL environment_end ( 'POST-PROC' )
276  !
277  CALL stop_pp()
278  !
279END PROGRAM pp
280