Projected band structures now works also in the

lsda case (in collaboration with B. Thakur).

src/manage_anhar.f90

@@ -28,7 +28,8 @@ SUBROUTINE manage_anhar()
 IMPLICIT NONE
 
 INTEGER :: itemp, igeom
-CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps
+CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps, &
+                      filepbs
 LOGICAL :: all_geometry_done, all_el_free, ldummy
 
 CALL check_all_geometries_done(all_geometry_done)
@@ -337,11 +338,11 @@ SUBROUTINE manage_anhar()
 !
 CALL write_gruneisen_band(flfrq_thermo,flvec_thermo)
 CALL set_files_for_plot(3, flfrq_thermo, filedata, filerap, &
-                                      fileout, gnu_filename, filenameps)
-CALL plotband_sub(3, filedata, filerap, fileout, gnu_filename, filenameps)
+                                      fileout, gnu_filename, filenameps, filepbs)
+CALL plotband_sub(3, filedata, filerap, fileout, gnu_filename, filenameps, filepbs)
 CALL set_files_for_plot(4, flfrq_thermo, filedata, filerap,  &
-                                       fileout, gnu_filename, filenameps)
-CALL plotband_sub(4, filedata, filerap, fileout, gnu_filename, filenameps)
+                                       fileout, gnu_filename, filenameps, filepbs)
+CALL plotband_sub(4, filedata, filerap, fileout, gnu_filename, filenameps, filepbs)
 !
 !    fit the frequencies of the dos mesh with a polynomial
 !
@@ -388,7 +389,8 @@ SUBROUTINE manage_anhar_anis()
 
 IMPLICIT NONE
 INTEGER :: itemp, itempp, igeom, startt, lastt, idata, ndata
-CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps
+CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps, &
+                      filepbs
 REAL(DP), ALLOCATABLE :: phf(:)
 LOGICAL :: all_geometry_done, all_el_free, ldummy
 INTEGER :: compute_nwork
@@ -680,8 +682,8 @@ SUBROUTINE manage_anhar_anis()
 
 CALL write_gruneisen_band_anis(flfrq_thermo,flvec_thermo)
 CALL set_files_for_plot(4, flfrq_thermo, filedata, filerap, &
-                                          fileout, gnu_filename, filenameps)
-CALL plotband_sub(4, filedata, filerap, fileout, gnu_filename, filenameps)
+                                          fileout, gnu_filename, filenameps, filepbs)
+CALL plotband_sub(4, filedata, filerap, fileout, gnu_filename, filenameps, filepbs)
 CALL plot_gruneisen_band_anis(flfrq_thermo)
 !
 !    fit the frequencies of the dos mesh with a polynomial

src/manage_bands.f90

@@ -27,7 +27,8 @@ SUBROUTINE manage_bands()
 IMPLICIT NONE
 
 INTEGER :: nspin0, exit_status, ierr
-CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps
+CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps, &
+                      filepbs
 CHARACTER(LEN=80)  :: message
 
 ierr=0
@@ -84,14 +85,14 @@ SUBROUTINE manage_bands()
    IF (is_a_path) THEN
       DO spin_component = 1, nspin0
          CALL set_files_for_plot(1, ' ', filedata, filerap, &
-                                           fileout, gnu_filename, filenameps)
+                                          fileout, gnu_filename, filenameps, filepbs)
          CALL plotband_sub(1,filedata, filerap, fileout, &
-                                           gnu_filename, filenameps)
+                                           gnu_filename, filenameps, filepbs)
       ENDDO
    ELSEIF (q2d) THEN
       spin_component=1
       CALL set_files_for_plot(1, ' ', filedata, filerap, &
-                                           fileout, gnu_filename, filenameps)
+                                           fileout, gnu_filename, filenameps, filepbs)
       CALL plot_ef(filedata, gnu_filename, filenameps)
    ENDIF
 ENDIF

src/manage_ph_postproc.f90

@@ -21,7 +21,8 @@ SUBROUTINE manage_ph_postproc(igeom)
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: igeom
 
-CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps
+CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps, &
+                      filepbs
 !
 !   Compute the interatomic force constants from the dynamical matrices
 !   written on file
@@ -33,9 +34,9 @@ SUBROUTINE manage_ph_postproc(igeom)
 IF (set_internal_path) CALL set_bz_path()
 CALL write_ph_dispersions()
 CALL set_files_for_plot(2, ' ', filedata, filerap, fileout, &
-                                                    gnu_filename, filenameps)
+                                            gnu_filename, filenameps, filepbs)
 IF (disp_nqs>0) CALL plotband_sub(2, filedata, filerap, fileout, &
-                                                    gnu_filename, filenameps)
+                                            gnu_filename, filenameps, filepbs)
 !
 !   Compute the harmonic thermodynamic quantities
 !

src/manage_surface_states.f90

@@ -27,12 +27,21 @@ SUBROUTINE manage_surface_states()
   INTEGER :: nks1, nks2, nks1tot, nks2tot
   LOGICAL :: exst
   INTEGER :: find_free_unit
-!
+  CHARACTER(LEN=256) :: fileprojlayer 
+  CHARACTER(LEN=6) :: int_to_char
+
+ !
 !   Find which k points must be done by this pool
 !
   CALL find_nks1nks2(1,nkstot,nks1tot,nks1,nks2tot,nks2,spin_component)
 
-  IF (ionode) INQUIRE( FILE = TRIM(flprojlayer), EXIST = exst )
+  fileprojlayer = flprojlayer
+
+  IF (nspin==2) &
+         fileprojlayer = TRIM(flprojlayer)// &
+                              '.'//TRIM(int_to_char(spin_component))
+
+  IF (ionode) INQUIRE( FILE = TRIM(fileprojlayer), EXIST = exst )
   CALL mp_bcast(exst, ionode_id, intra_image_comm)
 !
 !   the file with the projections is created here if it does not exist,
@@ -47,9 +56,9 @@ SUBROUTINE manage_surface_states()
                                                                 surface2)
   IF (ionode) THEN
      iun=find_free_unit()
-     OPEN(UNIT=iun, FILE=TRIM(flprojlayer), STATUS='unknown', ERR=400, &
+     OPEN(UNIT=iun, FILE=TRIM(fileprojlayer), STATUS='unknown', ERR=400, &
                                                             IOSTAT=ios)
-     WRITE(iun, '(5i8)') nat, nlayers, nbnd, nkstot, nspin     
+     WRITE(iun, '(5i8)') nat, nlayers, nbnd, nks2tot-nks1tot+1, nspin     
      WRITE(iun, '(4i8)') surface1, surface2    
      DO ik=nks1tot, nks2tot
         DO ibnd=1, nbnd

src/plot_gruneisen_band_anis.f90

@@ -28,7 +28,7 @@ SUBROUTINE plot_gruneisen_band_anis(flinput)
 CHARACTER(LEN=256), INTENT(IN) :: flinput
 CHARACTER(LEN=256) :: filename, save_flpsgrun, save_flgrun, save_flgnuplot, &
                       save_flpgrun
-CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps
+CHARACTER(LEN=256) :: filedata, filerap, fileout, gnu_filename, filenameps, filepbs
 
 CHARACTER(LEN=6), EXTERNAL :: int_to_char
 INTEGER :: nvar, icrys
@@ -48,8 +48,8 @@ SUBROUTINE plot_gruneisen_band_anis(flinput)
    flpgrun = TRIM(save_flpgrun)//'_'//TRIM(int_to_char(icrys))
    flgnuplot = TRIM(save_flgnuplot)//'_'//TRIM(int_to_char(icrys))
    CALL set_files_for_plot(3, flinput, filedata, filerap, &
-                                       fileout, gnu_filename, filenameps)
-   CALL plotband_sub(3,filedata, filerap, fileout, gnu_filename, filenameps)
+                                       fileout, gnu_filename, filenameps, filepbs)
+   CALL plotband_sub(3,filedata, filerap, fileout, gnu_filename, filenameps, filepbs)
 END DO
 
 flgrun=save_flgrun

src/plotband_sub.f90

@@ -9,7 +9,7 @@
 !
 !----------------------------------------------------------------------
 SUBROUTINE plotband_sub(icode, filedata, filerap, fileout, &
-                                      gnu_filename, filenameps )
+                                      gnu_filename, filenameps, filepbs)
 !----------------------------------------------------------------------
   !
   ! reads data files produced by "bands_sub", produces
@@ -51,7 +51,7 @@ SUBROUTINE plotband_sub(icode, filedata, filerap, fileout, &
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: icode
   CHARACTER(LEN=256), INTENT(IN) :: filedata, filerap, fileout, &
-                                               gnu_filename, filenameps
+                                         gnu_filename, filenameps, filepbs
 !
 !  path variables
 !
@@ -832,7 +832,7 @@ SUBROUTINE plotband_sub(icode, filedata, filerap, fileout, &
   IF (lprojpbs) CALL proj_band_structure(kx, e_eff, tot_points, nbnd, &
                     ymin, ymax, eref, e_rap, nrap, nbnd_rapk, start_rapk, &
                     nlines, start_point_eff, last_point_eff, &
-                    nrap_plot_eff, rap_plot_eff )
+                    nrap_plot_eff, rap_plot_eff, filepbs)
 
   IF (identify_sur) CALL plot_surface_states(nbnd, tot_points, nlines, kx, &
                            e_rap, ymin, ymax, eref, nrap, nbnd_rapk, &
@@ -901,7 +901,8 @@ SUBROUTINE plotband_sub(icode, filedata, filerap, fileout, &
   IF (identify_sur) THEN
      DEALLOCATE (lsurface_state_eff)
      DEALLOCATE (lsurface_state_rap)
-  ENDIF
+     DEALLOCATE (lsurface_state)
+   ENDIF
 
 !
 !  last deallocate the variables that define the path read from file

src/proj_band_structure.f90

@@ -9,7 +9,7 @@
 !------------------------------------------------------------
  SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
                  e_rap, nrap, nbnd_rapk, start_rapk, nlines_, start_point_, &
-                 last_point_, nrap_plot, rap_plot )
+                 last_point_, nrap_plot, rap_plot, filepbs)
 !------------------------------------------------------------
 !
 !     This routine uses gnuplot library to plot 
@@ -69,6 +69,7 @@ SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
 USE gnuplot,  ONLY : gnuplot_polygon, gnuplot_line, gnuplot_write_command
 USE point_group, ONLY : convert_rap
 USE noncollin_module, ONLY : lspinorb
+USE control_thermo,   ONLY : spin_component
 USE io_global, ONLY : ionode, ionode_id, stdout
 USE mp_images, ONLY : intra_image_comm
 USE mp, ONLY : mp_bcast
@@ -81,6 +82,7 @@ SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
                        nbnd_rapk(12,nks_), start_rapk(12,nks_)
 INTEGER, INTENT(IN) :: nrap_plot(nks_), rap_plot(12,nks_)
 REAL(DP), INTENT(IN) :: emin, emax, eref_
+CHARACTER(LEN=256) :: filepbs
 REAL(DP) :: e(nbnd_,nks_), kx(nks_), e_rap(nbnd_, nks_)
 
 INTEGER :: nks   ! this is the number of k point of a single path
@@ -180,7 +182,7 @@ SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
 
    IF (ionode) THEN
       iun=find_free_unit()
-      OPEN(UNIT=iun,FILE=TRIM(flpbs),STATUS='unknown',ERR=100,IOSTAT=ios)
+      OPEN(UNIT=iun,FILE=TRIM(filepbs),STATUS='unknown',ERR=100,IOSTAT=ios)
       WRITE(iun, '(3i5,f12.6)') nbnd, nks, nlines, eref
       DO ilines=1, nlines
          WRITE(iun,'(2i5)') start_point(ilines), last_point(ilines)
@@ -202,7 +204,7 @@ SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
 !  In this run the PBS information is read from file
 !
    IF (ionode) THEN
-      OPEN(UNIT=iun, FILE=TRIM(flpbs), STATUS='old', ERR=200, IOSTAT=ios)     
+      OPEN(UNIT=iun, FILE=TRIM(filepbs), STATUS='old', ERR=200, IOSTAT=ios)     
       READ(iun, '(3i5,f12.6)') nbnd, nks, nlines, eref
    ENDIF
 200 CALL mp_bcast(ios,ionode_id,intra_image_comm)
@@ -247,6 +249,7 @@ SUBROUTINE proj_band_structure(kx, e, nks_, nbnd_, emin, emax, eref_, &
 !     
 ! et1 and et2 (min/max/min/max...) is in a single variable eth(2*nbnd,nks)
 !
+
 ALLOCATE(eth(2*nbnd,nks))
 ALLOCATE(lth(2*nbnd,nks))
 ALLOCATE(x(nks))

src/read_state_densities.f90

@@ -14,6 +14,8 @@ SUBROUTINE read_state_densities()
 USE control_2d_bands,  ONLY : averag, vacuum, nlayers, identify_sur, &
                               surface1, surface2
 USE data_files,        ONLY : flprojlayer
+USE control_thermo,     ONLY : spin_component
+USE lsda_mod,           ONLY : nspin_=>nspin
 USE io_global,         ONLY : ionode, ionode_id
 USE mp_images,         ONLY : intra_image_comm
 USE mp,                ONLY : mp_bcast
@@ -24,16 +26,23 @@ SUBROUTINE read_state_densities()
 INTEGER :: iun, ios, idum, ilayer, ik, ibnd, nspin, nat_, nbnd_, nkstot_
 INTEGER :: ispin
 INTEGER :: find_free_unit
+CHARACTER(LEN=256) :: fileprojlayer
+CHARACTER(LEN=6) :: int_to_char
 
 IF (identify_sur) THEN
-   IF (ionode) &
-        INQUIRE( FILE = TRIM(flprojlayer), EXIST = exst ) 
+
+      fileprojlayer = flprojlayer
+      IF (nspin_==2) &
+         fileprojlayer = TRIM(flprojlayer)// &
+                              '.'//TRIM(int_to_char(spin_component))
+      IF (ionode) &
+        INQUIRE( FILE = TRIM(fileprojlayer), EXIST = exst ) 
      CALL mp_bcast(exst,ionode_id,intra_image_comm)
 
      IF (exst) THEN
         iun=find_free_unit()
         IF (ionode) THEN
-           OPEN(UNIT=iun,FILE=TRIM(flprojlayer),STATUS='old',ERR=300,&
+          OPEN(UNIT=iun,FILE=TRIM(fileprojlayer),STATUS='old',ERR=300,&
                                                               IOSTAT=ios)
            READ(iun, '(5i8)') nat, nlayers, nbnd_, nkstot_, nspin     
         ENDIF

src/set_files_names.f90

@@ -179,7 +179,7 @@ END SUBROUTINE initialize_el_file_names
 
 !----------------------------------------------------------------
 SUBROUTINE set_files_for_plot(icode, file_disp, filedata, filerap, fileout, &
-                                     gnu_filename, filenameps)
+                                     gnu_filename, filenameps, filepbs)
 !----------------------------------------------------------------
 !
 !   This routine receives as input a code of what we want to plot
@@ -202,18 +202,19 @@ SUBROUTINE set_files_for_plot(icode, file_disp, filedata, filerap, fileout, &
 USE control_thermo,   ONLY : spin_component
 USE lsda_mod,         ONLY : nspin
 
-USE data_files,       ONLY : flpgrun, flpband, filband, flfrq, flgrun
+USE data_files,       ONLY : flpgrun, flpband, filband, flfrq, flgrun, &
+                             flpbs
 USE postscript_files, ONLY : flpsband, flpsdisp, flpsgrun
 USE control_gnuplot,  ONLY : flgnuplot, flext
-
+USE control_2d_bands, ONLY : lprojpbs
 USE io_global,        ONLY : stdout
 
 IMPLICIT NONE
 
 INTEGER, INTENT(IN) :: icode
 CHARACTER(LEN=256), INTENT(IN) :: file_disp
 CHARACTER(LEN=256), INTENT(OUT) :: filedata, filerap, fileout, gnu_filename, &
-                                   filenameps 
+                                   filenameps, filepbs 
 CHARACTER(LEN=6) :: int_to_char
 !
 !  first the file with the data 
@@ -247,10 +248,14 @@ SUBROUTINE set_files_for_plot(icode, file_disp, filedata, filerap, fileout, &
         fileout=' '
      ELSEIF (icode==1) THEN
         fileout="band_files/"//TRIM(flpband)
-        IF (nspin==2) &
+        IF (nspin==2) THEN
            fileout="band_files/"//TRIM(flpband)//"."//&
-                                  TRIM(int_to_char(spin_component))
-     ELSEIF (icode==2) THEN
+                   TRIM(int_to_char(spin_component))
+           IF (lprojpbs) &
+           filepbs =TRIM(flpbs)//"."//&
+                    TRIM(int_to_char(spin_component))
+        ENDIF
+   ELSEIF (icode==2) THEN
         fileout="phdisp_files/"//TRIM(flpband)
      ENDIF
   ELSEIF (icode==3) THEN

src/surface.f90

@@ -241,6 +241,8 @@ SUBROUTINE identify_surface_states(nat, nbnd, nkstot, e, rap)
 
 DEALLOCATE(plot)
 DEALLOCATE(sumna)
+DEALLOCATE(averag)
+DEALLOCATE(vacuum)
 
 RETURN
 END SUBROUTINE identify_surface_states