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