diff --git a/src/postw90/comms.F90 b/src/postw90/comms.F90
index a9cd4975f..450d1c743 100644
--- a/src/postw90/comms.F90
+++ b/src/postw90/comms.F90
@@ -99,6 +99,7 @@ module w90_comms
   interface comms_gatherv
 !     module procedure comms_gatherv_int    ! to be done
      module procedure comms_gatherv_real
+     module procedure comms_gatherv_logical
 !     module procedure comms_gatherv_cmplx
   end interface comms_gatherv
 
@@ -846,6 +847,36 @@ subroutine comms_gatherv_real(array,localcount,rootglobalarray,counts,displs)
 
   end subroutine comms_gatherv_real
 
+  subroutine comms_gatherv_logical(array,localcount,rootglobalarray,counts,displs)
+    !! Gather real data to root node
+    implicit none
+
+    logical, intent(inout)           :: array
+    !! local array for sending data
+    integer, intent(in)                       :: localcount
+    !! localcount elements will be sent to the root node
+    logical, intent(inout)           :: rootglobalarray
+    !! array on the root node to which data will be sent
+    integer, dimension(num_nodes), intent(in) :: counts
+    !! how data should be partitioned, see MPI documentation or
+    !! function comms_array_split
+    integer, dimension(num_nodes), intent(in) :: displs
+
+#ifdef MPI
+    integer :: error
+
+    call MPI_gatherv(array,localcount,MPI_logical,rootglobalarray,counts,&
+         displs,MPI_logical,root_id,mpi_comm_world,error)
+
+    if(error.ne.MPI_success) then
+       call io_error('Error in comms_gatherv_logical')
+    end if
+#else
+!    rootglobalarray(1:localcount)=array(1:localcount)
+#endif
+
+  end subroutine comms_gatherv_logical
+
   subroutine comms_scatterv_real(array,localcount,rootglobalarray,counts,displs)
     !! Scatter data from root node
     implicit none
diff --git a/src/postw90/kslice.F90 b/src/postw90/kslice.F90
index 624e1fd31..fc44fb3ac 100644
--- a/src/postw90/kslice.F90
+++ b/src/postw90/kslice.F90
@@ -25,10 +25,6 @@ module w90_kslice
   !!
   !!    kslice_corner(1:3) is the lower left corner 
   !!    kslice_b1(1:3) and kslice_b2(1:3) are the vectors subtending the slice
-  !!
-  !!---------------------------------
-  !! TO DO: Parallelize over k-points
-  !!---------------------------------
 
   implicit none
 
@@ -62,9 +58,10 @@ subroutine k_slice
     use w90_berry, only          : berry_get_imf_klist,berry_get_imfgh_klist
     use w90_constants, only      : bohr
 
-    integer           ::  itot,i1,i2,n,n1,n2,n3,i
-    integer           :: zdataunit,coorddataunit,& 
-                         bandsunit,scriptunit,dataunit
+    integer, dimension(0:num_nodes-1) :: counts, displs
+
+    integer           :: iloc,itot,i1,i2,n,n1,n2,n3,i,nkpts,my_nkpts
+    integer           :: scriptunit
     real(kind=dp)     :: avec_2d(3,3),avec_3d(3,3),bvec(3,3),yvec(3),zvec(3),&
                          b1mod,b2mod,ymod,cosb1b2,kcorner_cart(3),&
                          areab1b2,cosyb2,kpt(3),kpt_x,kpt_y,k1,k2,&
@@ -82,23 +79,27 @@ subroutine k_slice
     complex(kind=dp), allocatable :: UU(:,:)
     real(kind=dp),    allocatable :: eig(:)
 
-    ! Everything is done on the root node. However, we still have to
-    ! read and distribute the data if we are in parallel, so calls to
-    ! get_oper are done on all nodes
-   
-    plot_fermi_lines=.false.
-    if(index(kslice_task,'fermi_lines')>0) plot_fermi_lines=.true.
-    plot_curv=.false.
-    if(index(kslice_task,'curv')>0) plot_curv=.true.
-    plot_morb=.false.
-    if(index(kslice_task,'morb')>0) plot_morb=.true.
-    fermi_lines_color=.false.
-    if(index(kslice_fermi_lines_colour,'spin')>0) fermi_lines_color=.true.
-    heatmap=.false.
-    if(plot_curv .or. plot_morb) heatmap=.true.
-    if(plot_fermi_lines .and. fermi_lines_color .and. heatmap)&
+    ! Output data buffers
+    real(kind=dp), allocatable    :: coords(:,:),    my_coords(:,:), &
+                                     spndata(:,:),   my_spndata(:,:), &
+                                     bandsdata(:,:), my_bandsdata(:,:), &
+                                     zdata(:,:),     my_zdata(:,:)
+    logical, allocatable          :: spnmask(:,:),   my_spnmask(:,:)
+
+    plot_fermi_lines  = index(kslice_task,'fermi_lines') > 0
+    plot_curv         = index(kslice_task,'curv') > 0
+    plot_morb         = index(kslice_task,'morb') > 0
+    fermi_lines_color = kslice_fermi_lines_colour /= 'none'
+    heatmap           = plot_curv .or. plot_morb
+    if(plot_fermi_lines .and. fermi_lines_color .and. heatmap) then
          call io_error('Error: spin-colored Fermi lines not allowed in '&
          //'curv/morb heatmap plots') 
+    end if
+
+    if(on_root) then
+       call kslice_print_info(plot_fermi_lines, fermi_lines_color, &
+                              plot_curv, plot_morb)
+    end if
 
     call get_HH_R
     if(plot_curv.or.plot_morb) call get_AA_R
@@ -108,8 +109,6 @@ subroutine k_slice
     endif
     if(fermi_lines_color) call get_SS_R
 
-    if(on_root) then
-
        ! Set Cartesian components of the vectors (b1,b2) spanning the slice
        !
        bvec(1,:)=matmul(kslice_b1(:),recip_lattice(:,:))
@@ -149,46 +148,14 @@ subroutine k_slice
          square='False'
        end if  
 
-       write(stdout,'(/,/,1x,a)')&
-            'Properties calculated in module  k s l i c e'
-       write(stdout,'(1x,a)')&
-            '--------------------------------------------'
-
-       if(plot_fermi_lines) then
-          if(nfermi/=1) call io_error(&
-               'Must specify one Fermi level when kslice_task=fermi_lines')
-          select case(fermi_lines_color)
-          case(.false.)
-             write(stdout,'(/,3x,a)') '* Fermi lines'
-          case(.true.)
-             write(stdout,'(/,3x,a)') '* Fermi lines coloured by spin'
-          end select
-          write(stdout,'(/,7x,a,f10.4,1x,a)')&
-               '(Fermi level: ',fermi_energy_list(1),'eV)'
-       endif
-       if(plot_curv) then
-          if(berry_curv_unit=='ang2') then
-             write(stdout,'(/,3x,a)') '* Negative Berry curvature in Ang^2'
-          elseif(berry_curv_unit=='bohr2') then
-             write(stdout,'(/,3x,a)') '* Negative Berry curvature in Bohr^2'
-          endif
-          if(nfermi/=1) call io_error(&
-               'Must specify one Fermi level when kslice_task=curv')
-       elseif(plot_morb) then
-          write(stdout,'(/,3x,a)')&
-               '* Orbital magnetization k-space integrand in eV.Ang^2'
-          if(nfermi/=1) call io_error(&
-               'Must specify one Fermi level when kslice_task=morb')
-       endif
+       nkpts = (kslice_2dkmesh(1)+1)*(kslice_2dkmesh(2)+1)
 
-       write(stdout,'(/,/,1x,a)') 'Output files:' 
+       ! Partition set of k-points into junks
+       call comms_array_split(nkpts, counts, displs);
+       my_nkpts = counts(my_node_id)
 
-       if(.not.fermi_lines_color) then
-          coorddataunit=io_file_unit() 
-          filename=trim(seedname)//'-kslice-coord.dat'
-          write(stdout,'(/,3x,a)') filename
-          open(coorddataunit,file=filename,form='formatted')
-       endif
+       allocate(my_coords(2,my_nkpts))
+       if(heatmap) allocate(my_zdata(3,my_nkpts))
 
        if(plot_fermi_lines) then
           allocate(HH(num_wann,num_wann))
@@ -196,47 +163,19 @@ subroutine k_slice
           allocate(eig(num_wann))
           if(fermi_lines_color) then
              allocate(delHH(num_wann,num_wann,3))
-             dataunit=io_file_unit()
-             filename=trim(seedname)//'-kslice-fermi-spn.dat'
-             write(stdout,'(/,3x,a)') filename
-             open(dataunit,file=filename,form='formatted')
+             allocate(my_spndata(num_wann,my_nkpts))
+             allocate(my_spnmask(num_wann,my_nkpts))
+             my_spnmask = .false.
           else
-             bandsunit=io_file_unit()
-             filename=trim(seedname)//'-kslice-bands.dat'
-             write(stdout,'(/,3x,a)') filename
-             open(bandsunit,file=filename,form='formatted')
-             if(.not.heatmap) then
-                allocate(bnddataunit(num_wann))
-                do n=1,num_wann
-                   n1=n/100
-                   n2=(n-n1*100)/10
-                   n3=n-n1*100-n2*10
-                   bnddataunit(n)=io_file_unit()
-                   filename=trim(seedname)//'-bnd_'&
-                        //achar(48+n1)//achar(48+n2)//achar(48+n3)//'.dat'
-                   write(stdout,'(/,3x,a)') filename
-                   open(bnddataunit(n),file=filename,form='formatted')
-                enddo
-             endif
+             allocate(my_bandsdata(num_wann,my_nkpts))
           endif
        endif
 
-       if(plot_curv) then
-          zdataunit=io_file_unit()
-          filename=trim(seedname)//'-kslice-curv.dat'
-          write(stdout,'(/,3x,a)') filename
-          open(zdataunit,file=filename,form='formatted')
-       elseif(plot_morb) then
-          zdataunit=io_file_unit()
-          filename=trim(seedname)//'-kslice-morb.dat'
-          write(stdout,'(/,3x,a)') filename
-          open(zdataunit,file=filename,form='formatted')
-       end if
-
-       ! Loop over uniform mesh of k-points covering the slice,
+       ! Loop over local portion of uniform mesh of k-points covering the slice,
        ! including all four borders
        !
-       do itot=0,(kslice_2dkmesh(1)+1)*(kslice_2dkmesh(2)+1)-1
+       do iloc=1,my_nkpts
+          itot=iloc-1+displs(my_node_id)
           i2=itot/(kslice_2dkmesh(1)+1) ! slow
           i1=itot-i2*(kslice_2dkmesh(1)+1) !fast
           ! k1 and k2 are the coefficients of the k-point in the basis
@@ -254,7 +193,8 @@ subroutine k_slice
           ! with x along x_vec=b1 and y along y_vec
           kpt_x=k1*b1mod+k2*b2mod*cosb1b2
           kpt_y=k2*b2mod*cosyb2
-          if(.not.fermi_lines_color) write(coorddataunit,'(2E16.8)') kpt_x,kpt_y
+
+          my_coords(:,iloc) = [kpt_x, kpt_y]
 
           if(plot_fermi_lines) then
              if(fermi_lines_color) then
@@ -272,27 +212,21 @@ subroutine k_slice
                 call pw90common_fourier_R_to_k(kpt,HH_R,HH,0)
                 call utility_diagonalize(HH,num_wann,eig,UU)
              endif
-             do n=1,num_wann
-                if(.not.fermi_lines_color) then
-                   ! For python
-                   write(bandsunit,'(E16.8)') eig(n)
-                   ! For gnuplot, using 'grid data' format
-                    if(.not.heatmap) then
-                       write(bnddataunit(n),'(3E16.8)') kpt_x,kpt_y,eig(n)
-                       if(i1==kslice_2dkmesh(1) .and. i2/=kslice_2dkmesh(2))&
-                            write (bnddataunit(n),*) ' '
-                    endif
-                 else
+
+             if(allocated(my_bandsdata)) then
+                my_bandsdata(:,iloc) = eig(:)
+             else
+                my_spndata(:,iloc) = spn_k(:)
+                do n=1,num_wann
                    ! vdum = dE/dk projected on the k-slice
                    zhat=zvec/sqrt(dot_product(zvec,zvec))
                    vdum(:)=del_eig(n,:)-dot_product(del_eig(n,:),zhat)*zhat(:)
                    Delta_E=sqrt(dot_product(vdum,vdum))*Delta_k
 !                   Delta_E=Delta_E*sqrt(2.0_dp) ! optimize this factor
-                   if(abs(eig(n)-fermi_energy_list(1))<Delta_E)&
-                        write(dataunit,'(3E16.8)') kpt_x,kpt_y,spn_k(n)
-                endif
-             enddo
-          endif
+                   my_spnmask(n,iloc) = abs(eig(n)-fermi_energy_list(1))<Delta_E
+                end do
+             end if
+          end if
 
           if(plot_curv) then
              call berry_get_imf_klist(kpt,imf_k_list)
@@ -301,10 +235,8 @@ subroutine k_slice
              curv(3)=sum(imf_k_list(:,3,1))
              if(berry_curv_unit=='bohr2') curv=curv/bohr**2   
              ! Print _minus_ the Berry curvature 
-             write(zdataunit,'(3E16.8)') -curv(:)
-          end if
-
-          if(plot_morb) then
+             my_zdata(:,iloc) = -curv(:)
+          else if(plot_morb) then
              call berry_get_imfgh_klist(kpt,imf_k_list,img_k_list,imh_k_list)
              Morb_k=img_k_list(:,:,1)+imh_k_list(:,:,1)&
                    -2.0_dp*fermi_energy_list(1)*imf_k_list(:,:,1)
@@ -312,34 +244,102 @@ subroutine k_slice
              morb(1)=sum(Morb_k(:,1))
              morb(2)=sum(Morb_k(:,2))
              morb(3)=sum(Morb_k(:,3))
-             write(zdataunit,'(3E16.8)') morb(:)
+             my_zdata(:,iloc) = morb(:)
           end if
 
-       end do !itot
-       
+       end do !iloc
+
+    ! Send results to root process
+    if(on_root) then
+       allocate(coords(2,nkpts))
+    else
+       allocate(coords(1,1))
+    end if
+    call comms_gatherv(my_coords(1,1), 2*my_nkpts, &
+                       coords(1,1), 2*counts, 2*displs)
+
+    if(allocated(my_spndata)) then
+       if(on_root) then
+          allocate(spndata(num_wann,nkpts))
+       else
+          allocate(spndata(1,1))
+       end if
+       call comms_gatherv(my_spndata(1,1), num_wann*my_nkpts, &
+                          spndata(1,1), num_wann*counts, num_wann*displs)
+    end if
+
+    if(allocated(my_spnmask)) then
+       if(on_root) then
+          allocate(spnmask(num_wann,nkpts))
+       else
+          allocate(spnmask(1,1))
+       end if
+       call comms_gatherv(my_spnmask(1,1), num_wann*my_nkpts, &
+                          spnmask(1,1), num_wann*counts, num_wann*displs)
+    end if
+
+    if(allocated(my_bandsdata)) then
+       if(on_root) then
+          allocate(bandsdata(num_wann,nkpts))
+       else
+          allocate(bandsdata(1,1))
+       end if
+       call comms_gatherv(my_bandsdata(1,1), num_wann*my_nkpts, &
+                          bandsdata(1,1), num_wann*counts, num_wann*displs)
+    end if
+
+    if(allocated(my_zdata)) then
+       if(on_root) then
+          allocate(zdata(3,nkpts))
+       else
+          allocate(zdata(1,1))
+       end if
+       call comms_gatherv(my_zdata(1,1), 3*my_nkpts, &
+                          zdata(1,1), 3*counts, 3*displs)
+    end if
+
+    ! Write output files
+    if(on_root) then
+       ! set kpt_x and kpt_y to last evaluated point
+       kpt_x = coords(1,nkpts)
+       kpt_y = coords(2,nkpts)
+
+       write(stdout,'(/,/,1x,a)') 'Output files:'
+
        if(.not.fermi_lines_color) then
-          write(coorddataunit,*) ' '
-          close(coorddataunit)
-       endif
+         filename=trim(seedname)//'-kslice-coord.dat'
+         call write_data_file(filename, '(2E16.8)', coords)
+       end if
        
-       if(heatmap) then
-          write(zdataunit,*) ' '
-          close(zdataunit)
-       endif
-       if(plot_fermi_lines) then
-          if(fermi_lines_color) then
-             close(dataunit)
-          else
-             write(bandsunit,*) ' '
-             close(bandsunit)
-             if(.not.heatmap) then
-                do n=1,num_wann
-                   write(bnddataunit(n),*) ' '
-                   close(bnddataunit(n))
-                enddo
-             endif
+       if(allocated(bandsdata)) then
+          ! For python
+          filename=trim(seedname)//'-kslice-bands.dat'
+          call write_data_file(filename, '(E16.8)', &
+                               reshape(bandsdata, [1, nkpts*num_wann]))
+
+          ! For gnuplot, using 'grid data' format
+          if(.not.heatmap) then
+             do n=1,num_wann
+                n1=n/100
+                n2=(n-n1*100)/10
+                n3=n-n1*100-n2*10
+                bnddataunit(n)=io_file_unit()
+                filename=trim(seedname)//'-bnd_'&
+                     //achar(48+n1)//achar(48+n2)//achar(48+n3)//'.dat'
+
+                call write_coords_file(filename, '(3E16.8)', coords, &
+                                       reshape(bandsdata(n,:), [1, 1, nkpts]), &
+                                       blocklen=kslice_2dkmesh(1)+1)
+             enddo
           endif
-       endif
+       end if
+
+       if(allocated(spndata)) then
+          filename=trim(seedname)//'-kslice-fermi-spn.dat'
+          call write_coords_file(filename, '(3E16.8)', coords, &
+                                 reshape(spndata, [1, num_wann, nkpts]), &
+                                 spnmask)
+       end if
 
        if(plot_fermi_lines .and. .not.fermi_lines_color .and. .not.heatmap) then
           !
@@ -607,6 +607,117 @@ subroutine k_slice
  
 end subroutine k_slice
 
+  !===========================================================!
+  !                   PRIVATE PROCEDURES
+  !===========================================================!
+
+  subroutine kslice_print_info(plot_fermi_lines, fermi_lines_color, plot_curv, plot_morb)
+    use w90_io,         only : stdout, io_error
+    use w90_parameters, only : nfermi, fermi_energy_list, berry_curv_unit
+
+    logical, intent(in)     :: plot_fermi_lines, fermi_lines_color, plot_curv, plot_morb
+
+    write(stdout,'(/,/,1x,a)')&
+         'Properties calculated in module  k s l i c e'
+    write(stdout,'(1x,a)')&
+         '--------------------------------------------'
+
+    if(plot_fermi_lines) then
+       if(nfermi/=1) call io_error(&
+            'Must specify one Fermi level when kslice_task=fermi_lines')
+       select case(fermi_lines_color)
+       case(.false.)
+          write(stdout,'(/,3x,a)') '* Fermi lines'
+       case(.true.)
+          write(stdout,'(/,3x,a)') '* Fermi lines coloured by spin'
+       end select
+       write(stdout,'(/,7x,a,f10.4,1x,a)')&
+            '(Fermi level: ',fermi_energy_list(1),'eV)'
+    endif
+    if(plot_curv) then
+       if(berry_curv_unit=='ang2') then
+          write(stdout,'(/,3x,a)') '* Negative Berry curvature in Ang^2'
+       elseif(berry_curv_unit=='bohr2') then
+          write(stdout,'(/,3x,a)') '* Negative Berry curvature in Bohr^2'
+       endif
+       if(nfermi/=1) call io_error(&
+            'Must specify one Fermi level when kslice_task=curv')
+    elseif(plot_morb) then
+       write(stdout,'(/,3x,a)')&
+            '* Orbital magnetization k-space integrand in eV.Ang^2'
+       if(nfermi/=1) call io_error(&
+            'Must specify one Fermi level when kslice_task=morb')
+    endif
+
+  end subroutine kslice_print_info
+
+  subroutine write_data_file(filename, fmt, data)
+     use w90_io,        only : io_error, stdout, io_file_unit
+     use w90_constants, only : dp
+
+     character(len=*), intent(in)  :: filename, fmt
+     real(kind=dp), intent(in)     :: data(:,:)
+
+     integer :: n, i, fileunit
+
+     write(stdout,'(/,3x,a)') filename
+     fileunit = io_file_unit()
+     open(fileunit,file=filename,form='formatted')
+
+     n = size(data,2)
+     do i=1,n
+        write(fileunit,fmt) data(:,i)
+     end do
+
+     write(fileunit,*) ''
+     close(fileunit)
+  end subroutine
+
+  subroutine write_coords_file(filename, fmt, coords, vals, mask, blocklen)
+     use w90_io,        only : io_error, stdout, io_file_unit
+     use w90_constants, only : dp
+
+     character(len=*), intent(in)  :: filename, fmt
+     real(kind=dp), intent(in)     :: coords(:,:), vals(:,:,:)
+     logical, intent(in), optional :: mask(:,:)
+     integer, intent(in), optional :: blocklen
+
+     integer :: n, m, i, j, fileunit, bl
+
+     write(stdout,'(/,3x,a)') filename
+     fileunit = io_file_unit()
+     open(fileunit,file=filename,form='formatted')
+
+     n = size(vals,3)
+     m = size(vals,2)
+
+     if(present(mask)) then
+        do i = 1,n
+           do j = 1,m
+              if(mask(j,i)) then
+                 write(fileunit, fmt) coords(:,i), vals(:,j,i)
+              end if
+           end do
+        end do
+     else
+        if(present(blocklen)) then
+           bl = blocklen
+        else
+           bl = n
+        end if
+
+        do i = 1,n
+           do j = 1,m
+              write(fileunit, fmt) coords(:,i), vals(:,j,i)
+           end do
+           if(mod(i,bl) == 0) then
+              write(fileunit, *) ''
+           end if
+        end do
+     end if
+     write(fileunit,*) ''
+     close(fileunit)
+  end subroutine
 
 subroutine script_common(scriptunit,areab1b2,square)
 
diff --git a/src/postw90/postw90_common.F90 b/src/postw90/postw90_common.F90
index e20911fdf..04c0c448e 100644
--- a/src/postw90/postw90_common.F90
+++ b/src/postw90/postw90_common.F90
@@ -298,7 +298,12 @@ subroutine pw90common_wanint_param_dist
     call comms_bcast(kpath_task,len(kpath_task)) 
     call comms_bcast(kpath_bands_colour,len(kpath_bands_colour)) 
     call comms_bcast(kslice,1) 
-    call comms_bcast(kslice_task,len(kslice_task)) 
+    call comms_bcast(kslice_task,len(kslice_task))
+    call comms_bcast(kslice_corner(1),3)
+    call comms_bcast(kslice_b1(1),3)
+    call comms_bcast(kslice_b2(1),3)
+    call comms_bcast(kslice_2dkmesh(1),2)
+    call comms_bcast(kslice_fermi_lines_colour, len(kslice_fermi_lines_colour))
     call comms_bcast(transl_inv,1) 
     call comms_bcast(num_elec_per_state,1)
     call comms_bcast(scissors_shift,1)