Merge lp:~nickpapior/siesta/4.1-ts into lp:siesta/4.1

[Nick Papior]

Enabled DM calculation for TBtrans. Also speeded up some of the Green function calculations in TBtrans.

Secondly, some bugs related to the pivoting schemes has been fixed. These fixes are for systems
where there is a gap between electrodes.

Lastly, user-defined energy points are now a reality in tbtrans.

[Nick Papior (nickpapior)]

review: Approve

=== modified file 'Docs/tbtrans.tex'
--- Docs/tbtrans.tex	2018-04-06 19:25:57 +0000
+++ Docs/tbtrans.tex	2018-05-04 19:08:45 +0000
@@ -124,7 +124,13 @@
   \item Projected transmission of eigenstates
 
   \item Orbital resolved ``bond-currents'' which may subsequently be
-  analyzed to yield actual bond-currents.
+  analyzed to yield actual bond-currents
+
+  \item Density matrices using the Green function and/or the
+  spectral density
+  
+  \item COOP and COHP curves using the Green function and/or the
+  spectral density.
   
 \end{itemize}
 
@@ -800,7 +806,30 @@
 \end{equation}
 where the sum is over all electrodes, $\mathbf G$ and $\mathbf A$ are
 the Green and spectral function, respectively. Note that typically
-$\rho_{\mathrm{bound-states}}=0$. 
+$\rho_{\mathrm{bound-states}}=0$.
+
+The below two options enables the calculation of the energy resolved
+density matrices. In effect they may be used to construct
+$\mathrm{LDOS}(E)$ profiles using \sisl.
+
+\begin{fdflogicalF}{TBT.DM!Gf}
+  \fdfdepend{TBT.Atoms!Device}
+
+  Calculate the energy and $k$-resolved density matrix for the Green
+  function. The density matrix may be used to construct 
+  real-space LDOS profiles.
+  
+\end{fdflogicalF}
+
+\begin{fdflogicalF}{TBT.DM!A}
+  \fdfdepend{TBT.Atoms!Device}
+
+  Calculate the energy and $k$-resolved density matrix for the
+  electrode spectral functions. The density matrix may be used to
+  construct real-space LDOS profiles.
+  
+\end{fdflogicalF}
+
 
 In addition to the DOS analysis of the Green and spectral functions,
 the Crystal Orbital Overlap Population and Crystal Orbital Hamilton
@@ -1326,9 +1355,10 @@
     Thus \emph{a} and \emph{b} are energies.
 
 
-    \option[points/delta]%
+    \option[points|delta|file]%
     \fdfindex{TBT.Contour!<>!points}%
     \fdfindex{TBT.Contour!<>!delta}%
+    \fdfindex{TBT.Contour!<>!file}%
 
     Define the number of integration points/energy separation.
     If specifying the number of points an integer should be supplied.
@@ -1336,6 +1366,31 @@
     If specifying the separation between consecutive points an energy
     should be supplied (e.g. \fdf*{0.01 eV}).
 
+    Optionally one may specify a file which contains the energy points
+    and their weights.
+
+    This file has the same formatting as the \sysfile{TBT.CC} output
+    with some optional inputs. Below is an example input file.
+    \begin{shellexample}
+# There are 2 different input options:
+# 1. Re[E] Im[E] W  (optional unit)
+# 2. Re[E] W        (optional unit) (imaginary part will be device Eta)
+# If the unit is specified on any line, all subsequent lines will use
+# the specified unit. Default unit is eV!
+# Empty lines and lines starting with # will be ignored.
+-0.5  0.1    # E = -0.5 eV, weight (for integrating current) of 0.1 eV
+-0.01 0.1 Ry # E = -0.01 Ry and weight 0.1 Ry
+-0.02 0.1    # E = -0.02 Ry (above unit continue) and weight 0.1 Ry
+-0.2  0.1 eV # E = -0.2 eV and weight 0.1 eV
+-0.2 1. 0.1  # E = -0.2 eV and 1. eV eta and weight 0.1 eV
+\end{shellexample}
+
+    If the file specified is \sysfile{TBT.CC} the same energy points
+    will be used. Note that the resulting \sysfile{TBT.nc} file does
+    not store the energies as complex numbers, thus one cannot
+    subsequently extract the $\eta$ value used for the individual
+    energy points.
+
     \option[method]%
     \fdfindex{TBT.Contour!<>!method}%
 
@@ -1360,6 +1415,9 @@
 
       \note has \fdf*{opt precision <>}.
 
+      \option[user]%
+      User defined input via a file.
+
     \end{fdfoptions}
 
     \option[opt]%
@@ -1377,14 +1435,28 @@
   TBT.Contours.Eta 0. eV
   %block TBT.Contours
     line
-  %endblock TBT.Contours
+  %endblock
   %block TBT.Contour.line
      from -2. eV to 2. eV
       delta 0.01 eV
         method mid-rule
-  %endblock TBT.Contour.line
+  %endblock
 \end{fdfexample}
 
+An example of input using a file (note that regular contour setups may
+be used together with file-inputs)
+\begin{fdfexample}
+  TBT.Contours.Eta 0. eV
+  %block TBT.Contours
+    file
+  %endblock
+  %block TBT.Contour.file
+     from 2. eV to 2.5 eV
+      file my_energies
+  %endblock
+\end{fdfexample}
+Note that the energy specifications are necessary (due to internal
+bookkeeping).
 
 \subsection{Chemical potentials}
 
@@ -1934,23 +2006,17 @@
 
   Specify the precision used for storing DOS in NetCDF4.
 
-  See \fdf{TBT.CDF!Precision}.
-  
 \end{fdfentry}
 
 \begin{fdfentry}{TBT.CDF!T.Precision}[string]<{<\fdf{TBT.CDF!Precision}>}>
 
   Specify the precision used for storing transmission function in NetCDF4.
 
-  See \fdf{TBT.CDF!Precision}.
-  
 \end{fdfentry}
 
 \begin{fdfentry}{TBT.CDF!T.Eig.Precision}[string]<{<\fdf{TBT.CDF!Precision}>}>
 
   Specify the precision used for storing transmission eigenvalues in NetCDF4.
-
-  See \fdf{TBT.CDF!Precision}.
   
 \end{fdfentry}
 
@@ -1958,7 +2024,23 @@
 
   Specify the precision used for storing orbital current in NetCDF4.
 
-  See \fdf{TBT.CDF!Precision}.
+  \note This is heavily advised to be in single precision as this may
+  easily use large amounts of disk-space if in double precision.
+  
+\end{fdfentry}
+
+\begin{fdfentry}{TBT.CDF!DM.Precision}[string]<{<\fdf{TBT.CDF!Precision}>}>
+
+  Specify the precision used for storing density matrices in NetCDF4.
+
+  \note This is heavily advised to be in single precision as this may
+  easily use large amounts of disk-space if in double precision.
+  
+\end{fdfentry}
+
+\begin{fdfentry}{TBT.CDF!COOP.Precision}[string]<{<\fdf{TBT.CDF!Precision}>}>
+
+  Specify the precision used for storing COOP and COHP curves in NetCDF4.
 
   \note This is heavily advised to be in single precision as this may
   easily use large amounts of disk-space if in double precision.
@@ -1978,7 +2060,9 @@
   NetCDF4 file using:
 \begin{shellexample}
    nccopy -d 3 siesta.TBT.nc newsiesta.TBT.nc
-\end{shellexample}
+ \end{shellexample}
+
+ \note one \emph{can not} do parallel I/O together with compression.
   
 \end{fdfentry}
 
@@ -1988,7 +2072,7 @@
   MPI processors this may increase performance. 
 
   \note this automatically sets the compression to 0 (one cannot
-  compress and perform parallel IO)
+  compress and perform parallel IO).
 
 \end{fdflogicalF}
 

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-04-24 09:14:21 +0000
+++ Src/Makefile	2018-05-04 19:08:45 +0000
@@ -944,13 +944,14 @@
 m_ts_contour_eq.o: m_gauss_fermi_17.o m_gauss_fermi_18.o m_gauss_fermi_19.o
 m_ts_contour_eq.o: m_gauss_fermi_20.o m_gauss_fermi_22.o m_gauss_fermi_24.o
 m_ts_contour_eq.o: m_gauss_fermi_26.o m_gauss_fermi_28.o m_gauss_fermi_30.o
-m_ts_contour_eq.o: m_gauss_fermi_inf.o m_gauss_quad.o m_integrate.o m_ts_aux.o
-m_ts_contour_eq.o: m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
+m_ts_contour_eq.o: m_gauss_fermi_inf.o m_gauss_quad.o m_integrate.o m_io.o
+m_ts_contour_eq.o: m_ts_aux.o m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
 m_ts_contour_eq.o: m_ts_io_contour.o m_ts_io_ctype.o parallel.o precision.o
 m_ts_contour_eq.o: units.o
-m_ts_contour_neq.o: m_gauss_quad.o m_integrate.o m_ts_aux.o m_ts_cctype.o
-m_ts_contour_neq.o: m_ts_chem_pot.o m_ts_electype.o m_ts_io_contour.o
-m_ts_contour_neq.o: m_ts_io_ctype.o parallel.o precision.o units.o
+m_ts_contour_neq.o: m_gauss_quad.o m_integrate.o m_io.o m_ts_aux.o
+m_ts_contour_neq.o: m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
+m_ts_contour_neq.o: m_ts_io_contour.o m_ts_io_ctype.o parallel.o precision.o
+m_ts_contour_neq.o: units.o
 m_ts_debug.o: class_Sparsity.o class_TriMat.o geom_helper.o parallel.o
 m_ts_debug.o: precision.o
 m_ts_dm_update.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
@@ -993,7 +994,7 @@
 m_ts_io.o: class_SpData2D.o class_Sparsity.o geom_helper.o m_io_s.o m_ncdf_io.o
 m_ts_io.o: m_os.o m_sparse.o parallel.o precision.o sys.o
 m_ts_io_contour.o: precision.o units.o
-m_ts_io_ctype.o: parallel.o precision.o units.o
+m_ts_io_ctype.o: m_io.o parallel.o precision.o units.o
 m_ts_iodm.o: class_OrbitalDistribution.o class_SpData2D.o class_Sparsity.o
 m_ts_iodm.o: m_io_s.o m_os.o parallel.o precision.o
 m_ts_kpoints.o: files.o find_kgrid.o kpoint_grid.o m_ts_global_vars.o

=== modified file 'Src/class_TriMat.T90'
--- Src/class_TriMat.T90	2018-03-26 08:06:50 +0000
+++ Src/class_TriMat.T90	2018-05-04 19:08:45 +0000
@@ -77,7 +77,7 @@
   public :: NEW_TYPE, print_type, init_val
 
   public :: val
-  public :: index, index_sub, part_index
+  public :: index, index_sub, index_block, part_index
   public :: cum_rows
   public :: nrows_g
   public :: parts, which_part

=== modified file 'Src/m_pivot.F90'
--- Src/m_pivot.F90	2017-11-13 07:41:00 +0000
+++ Src/m_pivot.F90	2018-05-04 19:08:45 +0000
@@ -30,10 +30,14 @@
   integer, parameter :: PVT_PCG = 7
   integer, parameter :: PVT_REV_PCG = 8
 #ifdef SIESTA__METIS
-  integer, parameter :: PVT_METIS = 9
+  integer, parameter :: PVT_METIS_NODEND = 9
 #endif
   integer, parameter :: PVT_CONNECT = 10
   integer, parameter :: PVT_REV_CONNECT = 11
+#ifdef SIESTA__METIS
+  integer, parameter :: PVT_METIS_PARTGRAPHKWAY = 12
+  integer, parameter :: PVT_METIS_PARTGRAPHRECURSIVE = 13
+#endif
 
 contains
 
@@ -120,9 +124,15 @@
        call rev_GGPS(n,n_nzs,ncol,l_ptr,l_col,lsub,pvt , priority = priority )
        pvt%name = 'rev-General-Gibbs-Poole-Stockmeyer'
 #ifdef SIESTA__METIS
-    else if ( method == PVT_METIS ) then
-       call metis_pvt(n,n_nzs,ncol,l_ptr,l_col,lsub,pvt, priority = priority)
-       pvt%name = 'metis'
+    else if ( method == PVT_METIS_NODEND ) then
+       call metis_NodeND_pvt(n,n_nzs,ncol,l_ptr,l_col,lsub,pvt, priority = priority)
+       pvt%name = 'metis-NodeND'
+    else if ( method == PVT_METIS_PARTGRAPHKWAY ) then
+       call metis_PartGraphKway_pvt(n,n_nzs,ncol,l_ptr,l_col,lsub,pvt, priority = priority)
+       pvt%name = 'metis-PartGraphKway'
+    else if ( method == PVT_METIS_PARTGRAPHRECURSIVE ) then
+       call metis_PartGraphRecursive_pvt(n,n_nzs,ncol,l_ptr,l_col,lsub,pvt, priority = priority)
+       pvt%name = 'metis-PartGraphRecursive'
 #endif
     else
        call die('m_pivot: Programming error, unknown method')

=== modified file 'Src/m_pivot_methods.F90'
--- Src/m_pivot_methods.F90	2018-04-07 19:13:42 +0000
+++ Src/m_pivot_methods.F90	2018-05-04 19:08:45 +0000
@@ -30,7 +30,9 @@
   public :: rev_connectivity_graph
 
 #ifdef SIESTA__METIS
-  public :: metis_pvt
+  public :: metis_PartGraphKway_pvt
+  public :: metis_PartGraphRecursive_pvt
+  public :: metis_NodeND_pvt
 #endif
 
   public :: bandwidth, profile
@@ -1446,7 +1448,87 @@
 
 
 #ifdef SIESTA__METIS
-  subroutine metis_pvt(n,nnzs,n_col,l_ptr,l_col,sub,pvt, priority)
+
+  subroutine metis_adjacency_graph(n,nnzs,n_col,l_ptr,l_col,sub, &
+      xadj, adjncy, w, priority)
+    use iso_c_binding, only: c_int
+    integer, intent(in) :: n, nnzs
+    integer, intent(in) :: n_col(n), l_ptr(n), l_col(nnzs)
+    type(tRgn), intent(in) :: sub
+    integer(c_int), intent(inout), allocatable :: xadj(:), adjncy(:), w(:)
+    integer, intent(in), optional :: priority(n)
+    
+    integer :: io, i, ptr, ind, nc, j, nadj
+
+    ! Allocate adjacency graphs
+    allocate(xadj(0:sub%n), w(sub%n) )
+    
+    ! First count adjacencies
+    xadj(0) = 0 ! + 1 : fortran-style
+    do i = 1 , sub%n
+
+      io = sub%r(i)
+      ptr = l_ptr(io)
+      nc = n_col(io)
+
+      ! Count number of elements in 
+      ! the sub-space
+      nadj = 0
+      do ind = ptr + 1 , ptr + nc
+        if ( in_rgn(sub,l_col(ind)) ) then
+          ! Skip "on-site" connections
+          if ( l_col(ind) /= io ) then
+            nadj = nadj + 1
+          end if
+        end if
+      end do
+
+      xadj(i) = xadj(i-1) + nadj
+      
+      if ( present(priority) ) then
+        w(i) = priority(io) + 1
+      else
+        w(i) = 1
+      end if
+      
+    end do
+
+    ! transfer to local adjacency graph
+    allocate( adjncy(xadj(sub%n)) )
+    
+    ! Create adjncy 
+    nadj = 0
+    do i = 1 , sub%n
+      
+      io = sub%r(i)
+      ptr = l_ptr(io)
+      nc = n_col(io)
+      
+      ! Count number of elements in 
+      ! the sub-space
+      do ind = ptr + 1 , ptr + nc
+        j = rgn_pivot(sub,l_col(ind))
+        if ( j > 0 ) then
+          ! Skip "on-site" connections
+          if ( l_col(ind) /= io ) then
+            nadj = nadj + 1
+            adjncy(nadj) = j - 1 ! + 1 : fortran style
+          end if
+        end if
+      end do
+      
+      if ( nadj /= xadj(i) ) then
+        print *,i,nadj, xadj(i)
+        call die('metis_adjacency_graph: Error in creating &
+            &adjacency graph.')
+      end if
+      
+    end do
+    
+  end subroutine metis_adjacency_graph
+
+  
+  subroutine metis_NodeND_pvt(n,nnzs,n_col,l_ptr,l_col,sub,pvt, priority)
     use iso_c_binding, only: c_int, c_ptr, c_loc
     integer, intent(in) :: n, nnzs
     integer, intent(in) :: n_col(n), l_ptr(n), l_col(nnzs)
@@ -1460,7 +1542,6 @@
     integer(c_int), allocatable, target :: w(:)
     integer(c_int) :: iret, nvtxs, opts(100) ! In 5.0.1 it is 40, but ...
     type(c_ptr) :: wp
-    integer :: nadj
     
     interface
        integer(c_int) function METIS_SetDefaultOptions(opts) &
@@ -1480,8 +1561,7 @@
        end function METIS_NodeND
     end interface
 
-    ! variables for the loop
-    integer :: io, i, ptr, ind, nc, j
+    integer :: i
 
     call rgn_delete(pvt)
 
@@ -1494,7 +1574,7 @@
     if ( iret /= 1 ) then
        opts(6) = 256 ! increase debug level and re-run for full dbg-lvl
        iret = METIS_setdefaultoptions(opts)
-       call die('metis_pvt: Error on initializing default options.')
+       call die('metis_NodeND: Error on initializing default options.')
     end if
 
     ! set options
@@ -1504,74 +1584,15 @@
     opts(8)  =  1 ! NCUT == Number of cuts
     opts(10) =  1 ! NO2HOP == does not use 2 hop
     opts(11) =  1 ! MINCONN == Explicitly minimize the maximum connectivity
-    opts(12) =  1 ! CONTIG == Forces contiguous 
+    opts(12) =  0 ! CONTIG == Forces contiguous 
     opts(13) =  1 ! COMPRESS == compress similar adjacency nodes
     opts(16) =  1 ! NSEPS(1) == tries in the separator
 
-    ! Allocate adjacency graphs
-    allocate(xadj(0:sub%n), perm(sub%n), iperm(sub%n) , w(sub%n) )
-
-    ! First count adjacencies
-    xadj(0) = 0 ! + 1 : fortran-style
-    do i = 1 , sub%n
-
-       io = sub%r(i)
-       ptr = l_ptr(io)
-       nc = n_col(io)
-
-       ! Count number of elements in 
-       ! the sub-space
-       nadj = 0
-       do ind = ptr + 1 , ptr + nc
-          if ( in_rgn(sub,l_col(ind)) ) then
-             ! Skip "on-site" connections
-             if ( l_col(ind) /= io ) then
-                nadj = nadj + 1
-             end if
-          end if
-       end do
-       
-       xadj(i) = xadj(i-1) + nadj
-
-       if ( present(priority) ) then
-          w(i) = priority(io) + 1
-       else
-          w(i) = 1
-       end if
-       
-    end do
-
-    ! transfer to local adjacency graph
-    allocate( adjncy(xadj(sub%n)) )
-
-    ! Create adjncy 
-    nadj = 0
-    do i = 1 , sub%n
-
-       io = sub%r(i)
-       ptr = l_ptr(io)
-       nc = n_col(io)
-
-       ! Count number of elements in 
-       ! the sub-space
-       do ind = ptr + 1 , ptr + nc
-          j = rgn_pivot(sub,l_col(ind))
-          if ( j > 0 ) then
-             ! Skip "on-site" connections
-             if ( l_col(ind) /= io ) then
-                nadj = nadj + 1
-                adjncy(nadj) = j - 1 ! + 1 : fortran style
-             end if
-          end if
-       end do
-       
-       if ( nadj /= xadj(i) ) then
-          print *,i,nadj, xadj(i)
-          call die('metis_pvt: Error in creating &
-               &adjacency graph.')
-       end if
-
-    end do
+    ! Setup the adjacency graph
+    call metis_adjacency_graph(n,nnzs,n_col,l_ptr,l_col,sub, &
+        xadj, adjncy, w, priority)
+    
+    allocate(perm(sub%n), iperm(sub%n))
 
     ! Call metis
     wp = c_loc(w(1))
@@ -1583,26 +1604,273 @@
        opts(6) = 256
        iret = METIS_NodeND(nvtxs, xadj, adjncy, wp, opts, perm, iperm)
        print *,iret
-       call die('pivot_method: metis, error in pivoting.')
+       call die('metis_NodeND: error in pivoting.')
     end if
 
     ! Clean-up
-    deallocate(xadj,adjncy,iperm)
+    deallocate(xadj,adjncy,w,iperm)
 
     call rgn_init(pvt,sub%n)
 
-
     ! Transfer pivoting to actual pivoting index
     do i = 1 , sub%n
        pvt%r(i) = sub%r(perm(i)+1) ! - 1 : fortran style
     end do
 
     ! Clean-up
-    deallocate(perm,w)
-
-    if ( pvt%n /= sub%n ) call die('metis: Error in algorithm')
-
-  end subroutine metis_pvt
+    deallocate(perm)
+
+    if ( pvt%n /= sub%n ) call die('metis_NodeND: Error in algorithm')
+
+  end subroutine metis_NodeND_pvt
+
+  subroutine metis_PartGraphKway_pvt(n,nnzs,n_col,l_ptr,l_col,sub,pvt, priority)
+    use iso_c_binding, only: c_int, c_ptr, c_loc, c_null_ptr
+    integer, intent(in) :: n, nnzs
+    integer, intent(in) :: n_col(n), l_ptr(n), l_col(nnzs)
+    type(tRgn), intent(in) :: sub
+    type(tRgn), intent(inout) :: pvt
+    integer, intent(in), optional :: priority(n)
+
+    ! METIS variables
+    integer(c_int), allocatable :: xadj(:), adjncy(:)
+    integer(c_int), allocatable :: part(:)
+    integer(c_int), allocatable, target :: w(:)
+    integer(c_int) :: iret, nvtxs, ncon, nparts, old_objval, objval, opts(100) ! In 5.0.1 it is 40, but ...
+    type(c_ptr) :: wp
+    
+    interface
+       integer(c_int) function METIS_SetDefaultOptions(opts) &
+            bind(C, name="METIS_SetDefaultOptions")
+         use iso_c_binding, only: c_int
+         implicit none
+         integer(c_int), dimension(*) :: opts
+       end function METIS_SetDefaultOptions
+       integer(c_int) function METIS_PartGraphKway(nvtxs,ncon,xadj,adjncy,vwgt, &
+           vsize,adjwgt,nparts,tpwgts,ubvec,opts,objval,part) bind(C, name="METIS_PartGraphKway")
+         use iso_c_binding, only: c_int, c_ptr
+         implicit none
+         integer(c_int) :: nvtxs, ncon, nparts, objval
+         integer(c_int), dimension(*) :: xadj, adjncy, opts, part
+         type(c_ptr), value :: vwgt, vsize, adjwgt, tpwgts, ubvec
+       end function METIS_PartGraphKway
+    end interface
+
+    ! variables for the loop
+    integer :: i, ip, j
+    integer :: old_bw, bw
+    type(tRgn) :: next_pvt
+
+    call rgn_delete(pvt)
+
+    ! The following does C-style indexing 
+    ! as the internal METIS structure is a simple offset
+
+!   call METIS_setdefaultoptions(opts)
+    iret = METIS_setdefaultoptions(opts)
+    ! METIS_OK == 1
+    if ( iret /= 1 ) then
+       opts(6) = 256 ! increase debug level and re-run for full dbg-lvl
+       iret = METIS_setdefaultoptions(opts)
+       call die('metis_PartGraphKway: Error on initializing default options.')
+    end if
+
+    ! set options
+    opts(2)  =  0 ! OBJTYPE == edge-cut minimization
+    opts(3)  =  1 ! CTYPE == Sorted heavy-edge matching
+!    opts(5)  =  1 ! RTYPE == Greedy-based cut and volume refinement
+    opts(7)  = 20 ! NITER(10) == Number of iterations
+    opts(8)  =  1 ! NCUTS == Number of cuts
+    opts(10) =  1 ! NO2HOP == does not use 2 hop
+    opts(11) =  1 ! MINCONN == Explicitly minimize the maximum connectivity
+    opts(12) =  0 ! CONTIG == Forces contiguous 
+
+    ! Allocate adjacency graphs
+    call metis_adjacency_graph(n,nnzs,n_col,l_ptr,l_col,sub, &
+        xadj, adjncy, w, priority)
+    
+    allocate(part(sub%n))
+
+    ! Initialize the pivoting array
+    call rgn_init(pvt,sub%n)
+    call rgn_init(next_pvt,sub%n)
+
+    iret = 1
+    nvtxs = sub%n
+    ncon = 1
+    nparts = 1
+    old_bw = huge(1)
+    wp = c_loc(w(1))
+    do nparts = 2 , min(nvtxs / 2 + 1, nvtxs)
+
+      ! Call metis
+      iret = METIS_PartGraphKway(nvtxs, ncon, xadj, adjncy, wp, &
+          C_NULL_PTR, C_NULL_PTR, &
+          nparts, & ! number of parts
+          C_NULL_PTR, C_NULL_PTR, opts, objval, part)
+
+      ! An error forces us to exit loop
+      if ( iret /= 1 ) exit
+
+      j = 0
+      do ip = 0, nparts - 1
+        do i = 1 , sub%n
+          if ( part(i) == ip ) then
+            j = j + 1
+            next_pvt%r(j) = sub%r(i)
+          end if
+        end do
+      end do
+
+      ! Transfer pivoting to actual pivoting index
+      bw = bandwidth(n,nnzs,n_col,l_ptr,l_col,next_pvt)
+      if ( bw == old_bw ) then
+        if ( profile(n,nnzs,n_col,l_ptr,l_col,next_pvt) &
+            < profile(n,nnzs,n_col,l_ptr,l_col,pvt) ) then
+          pvt%r(:) = next_pvt%r(:)
+        end if
+      else if ( bw < old_bw ) then
+        pvt%r(:) = next_pvt%r(:)
+      end if
+
+    end do
+
+    if ( iret /= 1 ) then ! METIS_OK == 1
+      print *,iret
+      call die('metis_PartGraphKway: error in pivoting.')
+    end if
+
+    ! Clean-up
+    deallocate(xadj,adjncy,w,part)
+    call rgn_delete(next_pvt)
+
+    if ( pvt%n /= sub%n ) call die('metis_PartGraphKway: Error in algorithm')
+
+  end subroutine metis_PartGraphKway_pvt
+
+  subroutine metis_PartGraphRecursive_pvt(n,nnzs,n_col,l_ptr,l_col,sub,pvt, priority)
+    use iso_c_binding, only: c_int, c_ptr, c_loc, c_null_ptr
+    integer, intent(in) :: n, nnzs
+    integer, intent(in) :: n_col(n), l_ptr(n), l_col(nnzs)
+    type(tRgn), intent(in) :: sub
+    type(tRgn), intent(inout) :: pvt
+    integer, intent(in), optional :: priority(n)
+
+    ! METIS variables
+    integer(c_int), allocatable :: xadj(:), adjncy(:)
+    integer(c_int), allocatable :: part(:)
+    integer(c_int), allocatable, target :: w(:)
+    integer(c_int) :: iret, nvtxs, ncon, nparts, old_objval, objval, opts(100) ! In 5.0.1 it is 40, but ...
+    type(c_ptr) :: wp
+    
+    interface
+       integer(c_int) function METIS_SetDefaultOptions(opts) &
+            bind(C, name="METIS_SetDefaultOptions")
+         use iso_c_binding, only: c_int
+         implicit none
+         integer(c_int), dimension(*) :: opts
+       end function METIS_SetDefaultOptions
+       integer(c_int) function METIS_PartGraphRecursive(nvtxs,ncon,xadj,adjncy,vwgt, &
+           vsize,adjwgt,nparts,tpwgts,ubvec,opts,objval,part) bind(C, name="METIS_PartGraphRecursive")
+         use iso_c_binding, only: c_int, c_ptr
+         implicit none
+         integer(c_int) :: nvtxs, ncon, nparts, objval
+         integer(c_int), dimension(*) :: xadj, adjncy, opts, part
+         type(c_ptr), value :: vwgt, vsize, adjwgt, tpwgts, ubvec
+       end function METIS_PartGraphRecursive
+    end interface
+
+    ! variables for the loop
+    integer :: i, ip, j
+    integer :: old_bw, bw
+    type(tRgn) :: next_pvt
+
+    call rgn_delete(pvt)
+
+    ! The following does C-style indexing 
+    ! as the internal METIS structure is a simple offset
+
+!   call METIS_setdefaultoptions(opts)
+    iret = METIS_setdefaultoptions(opts)
+    ! METIS_OK == 1
+    if ( iret /= 1 ) then
+       opts(6) = 256 ! increase debug level and re-run for full dbg-lvl
+       iret = METIS_setdefaultoptions(opts)
+       call die('metis_PartGraphRecursive: Error on initializing default options.')
+    end if
+
+    ! set options
+    opts(2)  =  0 ! OBJTYPE == edge-cut minimization
+    opts(3)  =  1 ! CTYPE == Sorted heavy-edge matching
+!    opts(5)  =  1 ! RTYPE == Greedy-based cut and volume refinement
+    opts(7)  = 20 ! NITER(10) == Number of iterations
+    opts(8)  =  1 ! NCUTS == Number of cuts
+    opts(10) =  1 ! NO2HOP == does not use 2 hop
+    opts(11) =  1 ! MINCONN == Explicitly minimize the maximum connectivity
+    opts(12) =  0 ! CONTIG == Forces contiguous 
+
+    ! Allocate adjacency graphs
+    call metis_adjacency_graph(n,nnzs,n_col,l_ptr,l_col,sub, &
+        xadj, adjncy, w, priority)
+    
+    allocate(part(sub%n))
+
+    ! Initialize the pivoting array
+    call rgn_init(pvt,sub%n)
+    call rgn_init(next_pvt,sub%n)
+
+    iret = 1
+    nvtxs = sub%n
+    ncon = 1
+    nparts = 1
+    old_bw = huge(1)
+    wp = c_loc(w(1))
+    do nparts = 2 , min(nvtxs / 2 + 1, nvtxs)
+
+      ! Call metis
+      iret = METIS_PartGraphRecursive(nvtxs, ncon, xadj, adjncy, wp, &
+          C_NULL_PTR, C_NULL_PTR, &
+          nparts, & ! number of parts
+          C_NULL_PTR, C_NULL_PTR, opts, objval, part)
+
+      ! An error forces us to exit loop
+      if ( iret /= 1 ) exit
+
+      j = 0
+      do ip = 0, nparts - 1
+        do i = 1 , sub%n
+          if ( part(i) == ip ) then
+            j = j + 1
+            next_pvt%r(j) = sub%r(i)
+          end if
+        end do
+      end do
+
+      ! Transfer pivoting to actual pivoting index
+      bw = bandwidth(n,nnzs,n_col,l_ptr,l_col,next_pvt)
+      if ( bw == old_bw ) then
+        if ( profile(n,nnzs,n_col,l_ptr,l_col,next_pvt) &
+            < profile(n,nnzs,n_col,l_ptr,l_col,pvt) ) then
+          pvt%r(:) = next_pvt%r(:)
+        end if
+      else if ( bw < old_bw ) then
+        pvt%r(:) = next_pvt%r(:)
+      end if
+
+    end do
+
+    if ( iret /= 1 ) then ! METIS_OK == 1
+      print *,iret
+      call die('metis_PartGraphRecursive: error in pivoting.')
+    end if
+
+    ! Clean-up
+    deallocate(xadj,adjncy,w,part)
+    call rgn_delete(next_pvt)
+
+    if ( pvt%n /= sub%n ) call die('metis_PartGraphRecursive: Error in algorithm')
+
+  end subroutine metis_PartGraphRecursive_pvt
 #endif
 
 
@@ -2126,30 +2394,27 @@
   function bandwidth(n,nnzs,n_col,l_ptr,l_col,sub) result(beta)
     integer, intent(in) :: n, nnzs, n_col(n), l_ptr(n), l_col(nnzs)
     type(tRgn), intent(in) :: sub
-    type(tRgn) :: s_sub, pvt
+    type(tRgn) :: pvt
     integer :: beta
     integer :: i, j, ind, idx
+
+    call rgn_init(pvt, n, val=0)
+    do i = 1, sub%n
+      pvt%r(sub%r(i)) = i
+    end do
+
     beta = 0
-
-    call rgn_copy(sub,s_sub)
-    call rgn_sort(s_sub)
-    call rgn_init(pvt,sub%n)
-    do i = 1 , sub%n
-       j = rgn_pivot(s_sub,sub%r(i))
-       pvt%r(j) = i
-    end do
-    
     do i = 1 , sub%n
        idx = sub%r(i)
        do ind = l_ptr(idx) + 1 , l_ptr(idx) + n_col(idx)
           ! figure out the pivoting place
-          j = rgn_pivot(s_sub,l_col(ind))
+          j = pvt%r(l_col(ind))
           if ( j <= 0 ) cycle
-          beta = max(beta,i-pvt%r(j))
+          beta = max(beta,i-j)
        end do
     end do
 
-    call rgn_delete(s_sub,pvt)
+    call rgn_delete(pvt)
 
   end function bandwidth
 
@@ -2157,32 +2422,29 @@
     integer, intent(in) :: n, nnzs, n_col(n), l_ptr(n), l_col(nnzs)
     type(tRgn), intent(in) :: sub
     integer(i8b) :: p
-    type(tRgn) :: s_sub, pvt
+    type(tRgn) :: pvt
     integer :: beta
     integer :: i, j, ind, idx
+
+    call rgn_init(pvt, n, val=0)
+    do i = 1, sub%n
+      pvt%r(sub%r(i)) = i
+    end do
+
     p = 0
-
-    call rgn_copy(sub,s_sub)
-    call rgn_sort(s_sub)
-    call rgn_init(pvt,sub%n)
-    do i = 1 , sub%n
-       j = rgn_pivot(s_sub,sub%r(i))
-       pvt%r(j) = i
-    end do
-    
     do i = 1 , sub%n
        idx = sub%r(i)
        beta = 0
        do ind = l_ptr(idx) + 1 , l_ptr(idx) + n_col(idx)
           ! figure out the pivoting place
-          j = rgn_pivot(s_sub,l_col(ind))
+          j = pvt%r(l_col(ind))
           if ( j <= 0 ) cycle
-          beta = max(beta,i-pvt%r(j))
+          beta = max(beta,i-j)
        end do
        p = p + beta
     end do
 
-    call rgn_delete(s_sub,pvt)
+    call rgn_delete(pvt)
 
   end function profile
 
@@ -2558,8 +2820,19 @@
           ! We will not follow the graph
           if ( lonly_sub ) exit
 
-          idx = idx_degree(D_LOW,n,nnzs,n_col,l_ptr,l_col,sub, skip = skip, &
+          ! Fake the skipping region
+          ! This is just to skip the first couple of zeros
+          call rgn_delete(con)
+          do idx = 1, skip%n
+            if ( skip%r(idx) > 0 ) then
+              con%n = skip%n - idx + 1
+              con%r => skip%r(idx:)
+              exit
+            end if
+          end do
+          idx = idx_degree(D_LOW,n,nnzs,n_col,l_ptr,l_col,sub, skip = con, &
                priority = priority)
+          call rgn_nullify(con)
           etr = sub%r(idx)
           
           ! Push the queue and the skip table
@@ -2576,7 +2849,7 @@
 
        ! Since the sort-degree is sorting into the same array, we have
        ! to ensure con to have the correct size
-       call rgn_init(con, sub%n - pvt%n)
+       call rgn_grow(con, sub%n - pvt%n)
 
        ! 3. Create the connectivity graph from idx (this will remove "back" 
        !    connected entries, hence no dublicates needs to be taken into 

=== modified file 'Src/m_ts_cctype.f90'
--- Src/m_ts_cctype.f90	2016-03-01 14:24:48 +0000
+++ Src/m_ts_cctype.f90	2018-05-04 19:08:45 +0000
@@ -61,6 +61,7 @@
   integer, parameter :: CC_BOOLE_MIX        = 103
   integer, parameter :: CC_MID              = 104
   integer, parameter :: CC_CONTINUED_FRAC   = 105
+  integer, parameter :: CC_USER             = 106
 
 
   ! Converts a method to a string format
@@ -116,7 +117,9 @@
     case ( CC_MID )
        str = 'Mid-rule'
     case ( CC_CONTINUED_FRAC )
-       str = 'Continued fraction'
+       str = 'Continued-fraction'
+    case ( CC_USER )
+       str = 'User'
     case default
        call die('Unknown method for the contour')
     end select
@@ -139,7 +142,9 @@
     case ( CC_MID )
        str = 'Mid-rule'
     case ( CC_CONTINUED_FRAC )
-       str = 'Continued fraction'
+       str = 'Continued-fraction'
+    case ( CC_USER )
+       str = 'User defined'
     case default
        call die('Unknown method for the contour')
     end select
@@ -175,6 +180,8 @@
        str = 'Mid-rule'
     case ( CC_CONTINUED_FRAC )
        str = 'Continued fraction'
+    case ( CC_USER )
+       str = 'User defined'
     case default
        call die('Unknown method for the contour')
     end select
@@ -203,6 +210,9 @@
          leqi(str,'continued-fraction') .or. &
          leqi(str,'cont-frac') ) then
        method = CC_CONTINUED_FRAC
+    else if ( leqi(str,'file') .or. &
+         leqi(str,'user') ) then
+       method = CC_USER
     else if ( leqi(str,'g-fermi') ) then
        method = CC_G_NF_0kT
        do i = G_NF_MIN_kT , G_NF_MAX_kT

=== modified file 'Src/m_ts_contour_eq.f90'
--- Src/m_ts_contour_eq.f90	2017-10-23 09:18:56 +0000
+++ Src/m_ts_contour_eq.f90	2018-05-04 19:08:45 +0000
@@ -442,7 +442,11 @@
        ! and not in "random" order
        idx = get_c_io_index(mu%Eq_seg(i))
 
-       if ( leqi(Eq_c(idx)%c_io%part,'circle') ) then
+       if ( leqi(Eq_c(idx)%c_io%part,'user') ) then
+
+          call contour_file(Eq_c(idx),mu,lift)
+
+        else if ( leqi(Eq_c(idx)%c_io%part,'circle') ) then
 
           call contour_Circle(Eq_c(idx),mu,R,cR)
 
@@ -513,7 +517,7 @@
 
       ! the radius can be calculated using two triangles in the circle
       ! there is no need to use the cosine-relations
-      R = .5_dp * cR / cos(alpha) ** 2
+      R = 0.5_dp * cR / cos(alpha) ** 2
 
       ! the real-axis center
       cR = a + R
@@ -1127,34 +1131,47 @@
           call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp)
 
        end if
-       
-       
+
+    case ( CC_USER )
+
+      ! Read the file information
+      call contour_file(c,mu,Eta)
+
     case default
        write(*,*) 'Method for contour ',trim(c%c_io%name), &
             ' could not be deciphered: ', c%c_io%method
        call die('Could not determine the line-integral')
     end select
-
-    ! I know this is "bad" practice, however, zero is a well-defined quantity.
-    set_c = sum(abs(c%c(:))) == 0._dp
-
+    
     ! get the index in the ID array (same index in w-array)
     call ID2idx(c,mu%ID,idx)
-
-    do i = 1 , c%c_io%N
-       if ( set_c ) then
+    
+    if ( method(c%c_io) == CC_USER ) then
+
+      do i = 1 , c%c_io%N
+        c%w(idx,i) = c%w(idx,i) * nf((real(c%c(i), dp) - mu%mu) / mu%kT)
+      end do
+      
+    else
+      
+      ! I know this is "bad" practice, however, zero is a well-defined quantity.
+      set_c = sum(abs(c%c(:))) == 0._dp
+
+      do i = 1 , c%c_io%N
+        if ( set_c ) then
           c%c(i) = dcmplx(ce(i),Eta)
-       else
+        else
           if ( abs(c%c(i) - dcmplx(ce(i),Eta)) > 1.e-10_dp ) then
-             call die('contour_line: Error on contour match')
+            call die('contour_line: Error on contour match')
           end if
-       end if
-
-       !ztmp = (c%c(i) - mu%mu) / mu%kT
-       c%w(idx,i) = cw(i) * nf((ce(i) - mu%mu) / mu%kT)
-
-    end do
-
+        end if
+
+        c%w(idx,i) = cw(i) * nf((ce(i) - mu%mu) / mu%kT)
+
+      end do
+      
+    end if
+    
     deallocate(ce,cw)
     
   end subroutine contour_line
@@ -1256,7 +1273,9 @@
 
     case default
 
-       ! we revert so that we can actually use the line-integral
+      ! we revert so that we can actually use the line-integral
+      ! The tail and line are equivalent in the sense that the
+      ! fermi functions are applied to the weights
        c%c_io%part = 'line'
 
        call contour_line(c,mu,Eta)
@@ -1438,6 +1457,110 @@
     
   end subroutine contour_continued_fraction
 
+  ! This routine will read the contour points from a file
+  subroutine contour_file(c,mu,Eta)
+    use m_io, only: io_assign, io_close
+    use fdf, only: fdf_convfac
+
+    type(ts_cw), intent(inout) :: c
+    type(ts_mu), intent(in) :: mu
+    ! The lifting into the complex plane
+    real(dp), intent(in) :: Eta
+
+    integer :: iu, iostat, ne, idx
+    logical :: exist
+    complex(dp) :: E , W
+    real(dp) :: rE, iE, rW, iW, conv
+    character(len=512) :: file, line
+    character(len=16) :: unit
+    
+    ! The contour type contains the file name in:
+    !  c%c_io%cN (weirdly enough)
+    file = c%c_io%cN
+    call ID2idx(c,mu%ID,idx)
+
+    call io_assign(iu)
+    
+    ! Open the contour file
+    inquire(file=trim(file), exist=exist)
+    if ( .not. exist ) then
+      call die('The file: '//trim(file)//' could not be found to read contour points!')
+    end if
+    
+    ! Open the file
+    open(iu, file=trim(file), form='formatted', status='old')
+    
+    ne = 0
+    ! The default unit is eV.
+    ! On every line an optional unit-specificer may be used to specify the
+    ! subsequent lines units (until another unit is specified)
+    conv = fdf_convfac('eV', 'Ry')
+    do
+      ! Now we have the line
+      read(iu, '(a)', iostat=iostat) line
+      if ( iostat /= 0 ) exit
+      if ( len_trim(line) == 0 ) cycle
+      line = trim(adjustl(line))
+      if ( line(1:1) == '#' ) cycle
+      
+      ! We have a line with energy and weight
+      ne = ne + 1
+      ! There are three optional ways of reading this
+      ! 1.  ReE ImE, ReW ImW [unit]
+      read(line, *, iostat=iostat) rE, iE, rW, iW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW, iW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(idx,ne) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 2.  ReE ImE, ReW [unit]
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, iE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(idx,ne) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 3.  ReE , ReW [unit]
+      iE = Eta
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE * conv,iE)
+        c%w(idx,ne) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      call die('Contour file: '//trim(file)//' is not formatted correctly. &
+          &Please read the documentation!')
+
+    end do
+
+    call io_close(iu)
+
+    if ( c%c_io%N /= ne ) then
+      call die('Error in reading the contour points from file: '//trim(file))
+    end if
+    
+  end subroutine contour_file
+  
   function Eq_E(id,step) result(c)
     integer, intent(in) :: id
     integer, intent(in), optional :: step

=== modified file 'Src/m_ts_contour_neq.f90'
--- Src/m_ts_contour_neq.f90	2017-08-22 15:09:30 +0000
+++ Src/m_ts_contour_neq.f90	2018-05-04 19:08:45 +0000
@@ -449,22 +449,26 @@
     type(ts_cw), intent(inout) :: c
     real(dp), intent(in) :: Eta
 
-    if ( leqi(c%c_io%part,'line') ) then
-       
-       call contour_line(c,Eta)
-       
+    if ( leqi(c%c_io%part,'user') ) then
+       
+      call contour_file(c,Eta)
+
+    else if ( leqi(c%c_io%part,'line') ) then
+
+      call contour_line(c,Eta)
+
     else if ( leqi(c%c_io%part,'tail') ) then
 
-       c%c_io%part = 'line'
-       
-       call contour_line(c,Eta)
-
-       c%c_io%part = 'tail'
+      c%c_io%part = 'line'
+
+      call contour_line(c,Eta)
+
+      c%c_io%part = 'tail'
 
     else
-       
-       call neq_die('Unrecognized contour type for the &
-            &non-equilibrium part.')
+
+      call neq_die('Unrecognized contour type for the &
+          &non-equilibrium part.')
        
     end if
     
@@ -602,6 +606,14 @@
        end if
 
        call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp)
+
+     case ( CC_USER )
+
+       call contour_file(c, Eta)
+
+       ! Immediately return as the user has specified *everything*
+       deallocate(ce, cw)
+       return
        
     case default
 
@@ -616,6 +628,109 @@
     
   end subroutine contour_line
 
+
+  ! This routine will read the contour points from a file
+  subroutine contour_file(c,Eta)
+    use m_io, only: io_assign, io_close
+    use fdf, only: fdf_convfac
+
+    type(ts_cw), intent(inout) :: c
+    ! The lifting into the complex plane
+    real(dp), intent(in) :: Eta
+
+    integer :: iu, iostat, ne
+    logical :: exist
+    complex(dp) :: E , W
+    real(dp) :: rE, iE, rW, iW, conv
+    character(len=512) :: file, line
+    character(len=16) :: unit
+    
+    ! The contour type contains the file name in:
+    !  c%c_io%cN (weirdly enough)
+    file = c%c_io%cN
+
+    call io_assign(iu)
+    
+    ! Open the contour file
+    inquire(file=trim(file), exist=exist)
+    if ( .not. exist ) then
+      call die('The file: '//trim(file)//' could not be found to read contour points!')
+    end if
+    
+    ! Open the file
+    open(iu, file=trim(file), form='formatted', status='old')
+    
+    ne = 0
+    ! The default unit is eV.
+    ! On every line an optional unit-specificer may be used to specify the
+    ! subsequent lines units (until another unit is specified)
+    conv = fdf_convfac('eV', 'Ry')
+    do
+      ! Now we have the line
+      read(iu, '(a)', iostat=iostat) line
+      if ( iostat /= 0 ) exit
+      if ( len_trim(line) == 0 ) cycle
+      line = trim(adjustl(line))
+      if ( line(1:1) == '#' ) cycle
+      
+      ! We have a line with energy and weight
+      ne = ne + 1
+      ! There are three optional ways of reading this
+      ! 1.  ReE ImE, ReW ImW [unit]
+      read(line, *, iostat=iostat) rE, iE, rW, iW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW, iW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 2.  ReE ImE, ReW [unit]
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, iE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 3.  ReE , ReW [unit]
+      iE = Eta
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE * conv,iE)
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      call die('Contour file: '//trim(file)//' is not formatted correctly. &
+          &Please read the documentation!')
+
+    end do
+
+    call io_close(iu)
+
+    if ( c%c_io%N /= ne ) then
+      call die('Error in reading the contour points from file: '//trim(file))
+    end if
+    
+  end subroutine contour_file
+
   function nEq_E(id,step) result(c)
     integer, intent(in) :: id
     integer, intent(in), optional :: step

=== modified file 'Src/m_ts_io_ctype.f90'
--- Src/m_ts_io_ctype.f90	2016-04-27 12:13:56 +0000
+++ Src/m_ts_io_ctype.f90	2018-05-04 19:08:45 +0000
@@ -264,6 +264,7 @@
 
     character(len=c_N) :: opt, val
     integer :: iS, iE
+    logical :: is_file
     
     ! if the block does not exist, return
     if ( .not. ts_exists_contour_block(prefix,suffix,bName,bfdf=bfdf) ) then
@@ -277,7 +278,7 @@
     c%name = trim(bName)
     
     ! We must ensure the block be organized as this:
-    !  part circle|line|tail
+    !  part circle|line|tail|user
     !  from <a> to <b> 
     !    points <N> / separation <d>
     !      method <method>
@@ -312,8 +313,10 @@
        c%part = 'line'
     else if ( leqi(c%part,'tail') ) then
        c%part = 'tail'
+    else if ( leqi(c%part,'user') ) then
+       c%part = 'user'
     else
-       call die('Part of the contour could not be recognized as circle|square|line|tail')
+       call die('Part of the contour could not be recognized as circle|square|line|tail|user')
     end if
    
 
@@ -361,17 +364,30 @@
     end if
     
     ! we now read the points or separation
+    is_file = .false.
     iS = search_fun('points',pline)
     if ( iS < 0 ) iS = search_fun('p',pline)
     iE = search_fun('delta',pline)
-    if ( iE < 0 ) iE = search_fun('sep',pline)
+    if ( iE < 0 ) iE = search_fun('d',pline)
     if ( iS < 0 .and. iE < 0 ) then
-       call die('Block: '//trim(bName)//' is not build correctly. &
-            &Could not decipher points/delta/separation')
+      iS = search_fun('file',pline)
+      if ( iS < 0 ) iS = search_fun('user',pline)
+      if ( iS < 0 ) then
+        call die('Block: '//trim(bName)//' is not build correctly. &
+            &Could not decipher points/delta/file')
+      end if
+      is_file = .true.
     end if
     
     ! if we have points we simply read in the number
-    if ( 0 <= iS ) then
+    if ( is_file ) then
+      
+      ! The file contains all the relevant information.
+      c%cN = fdf_bnames(pline, 2) ! file name
+      c%method = 'user' ! user-defined method
+      c%N = file_energy_points(c%cN) ! read number of energy points
+
+    else if ( 0 <= iS ) then
        c%N  = fdf_bintegers(pline,1,after=iS) ! first integer
        c%cN = characters(pline,1,-1,after=iS)
        if ( c%N < 1 ) then
@@ -397,17 +413,30 @@
     ! } "points"
 
     ! { "method <method>"
-    if ( .not. move2names() ) then
-       call die('Block: '//trim(bName)//'. &
+    if ( move2names() ) then
+      val = fdf_bnames(pline, 1)
+      if ( leqi(val, 'method') ) then
+        if ( fdf_bnnames(pline) < 2 ) then
+          call die('Block: '//trim(bName)//' has not described the method properly. &
+              &Must have method <method>')
+        end if
+
+        ! the method should be a one-name thing
+        c%method = fdf_bnames(pline,2)
+        
+      else if ( .not. is_file ) then
+        ! The method segment has not been found.
+        ! Error out
+        call die('Block: '//trim(bName)//'. &
             &Could not find method <method> segment in contour')
-    end if
-    if ( fdf_bnnames(pline) < 2 ) then
-       call die('Block: '//trim(bName)//' has not described the method properly. &
-            &Must have method <method>')
-    end if
-
-    ! the method should be a one-name thing
-    c%method = fdf_bnames(pline,2)
+      else if ( is_file ) then
+        ! The name is most probably an option, step back in the block
+        if ( .not. fdf_bbackspace(bfdf) ) then
+          call die('Block: '//trim(bName)//' parsing went wrong!')
+        end if
+      end if
+    end if
+    if ( is_file ) c%method = 'user'
 
     ! } "method"
 
@@ -432,7 +461,43 @@
     end do
     ! } "opt"
     
-  contains 
+  contains
+
+    function file_energy_points(file) result(ne)
+      use m_io, only: io_assign, io_close
+      
+      character(len=*), intent(in) :: file
+      integer :: iu
+      character(len=256) :: line
+      integer :: ne, iostat
+      logical :: exist
+
+      call io_assign(iu)
+
+      ! Open the contour file
+      inquire(file=trim(file), exist=exist)
+      if ( .not. exist ) then
+        call die('Block: '//trim(bName)//' requested an external contour file. &
+            &The file: '//trim(file)//' could not be found!')
+      end if
+
+      ! Open the file
+      open(iu, file=trim(file), form='formatted', status='old')
+
+      ne = 0
+      do
+        read(iu, '(a)', iostat=iostat) line
+        if ( iostat /= 0 ) exit
+        if ( len_trim(line) == 0 ) cycle
+        line = trim(adjustl(line))
+        if ( line(1:1) == '#' ) cycle
+        
+        ne = ne + 1
+      end do
+
+      call io_close(iu)
+      
+    end function file_energy_points
 
     function move2names() result(found)
       logical :: found
@@ -804,13 +869,17 @@
       
     write(*,'(4a)') '   from ',trim(c%ca),' to ',trim(c%cb)
 
-    if ( len_trim(c%cd) /= 0 ) then
-       ! we have delta designation
-       write(*,'(t7,a,tr1,a)') 'delta', trim(c%cd)
-
+    ! Ensure we correctly get the user-definition
+    if ( c%method == 'user' ) then
+      write(*,'(t7,a,tr1,a)') 'file', trim(c%cN)
+      
+    else if ( len_trim(c%cd) /= 0 ) then
+      ! we have delta designation
+      write(*,'(t7,a,tr1,a)') 'delta', trim(c%cd)
+      
     else       
-       ! Print the number of points...
-       write(*,'(t7,a,tr1,i0)') 'points', c%N
+      ! Print the number of points...
+      write(*,'(t7,a,tr1,i0)') 'points', c%N
       
     end if
 
@@ -831,6 +900,7 @@
     
   end subroutine ts_print_contour_block
 
+
   ! *****
   ! This routine fixes the inputs for the contours according to those given by 
   ! the input electrode

=== modified file 'Src/m_ts_method.f90'
--- Src/m_ts_method.f90	2018-02-27 21:03:50 +0000
+++ Src/m_ts_method.f90	2018-05-04 19:08:45 +0000
@@ -210,6 +210,10 @@
   subroutine set_type(typ,ia,na_u,lasto)
     integer, intent(in) :: typ, ia, na_u,lasto(0:na_u)
     integer :: i, no
+    if ( ia > na_u ) then
+      call die('Error in specifying the type of an atom!. &
+          &Atoms specified is above the total number of atoms!')
+    end if
     if ( a_type(ia) /= TYP_DEVICE ) then
       write(*,'(2(a,i0))') 'Trying to set atom ',ia,' to type: ',typ
       write(*,'(2(a,i0))') 'Atom ',ia,' is already: ',a_type(ia)

=== modified file 'Src/m_ts_options.F90'
--- Src/m_ts_options.F90	2018-04-06 19:25:57 +0000
+++ Src/m_ts_options.F90	2018-05-04 19:08:45 +0000
@@ -182,7 +182,7 @@
        ts_method = TS_MUMPS
 #endif
     else
-       call die('Unrecognized Transiesta solution method: '//trim(chars))
+       call die('Unrecognized TranSiesta solution method: '//trim(chars))
     end if
 
     ! currently this does not work
@@ -251,7 +251,7 @@
 
     ts_kT = fdf_get('TS.ElectronicTemperature',kT,'Ry')
     if ( ts_kT / Kelvin < 10._dp ) then
-       call die('transiesta electronic temperature *must* &
+       call die('TranSiesta electronic temperature *must* &
             &be larger than 10 kT')
     end if
     
@@ -277,7 +277,7 @@
        ! We do not allow the electronic temperature
        ! to be below 10 kT
        if ( mus(i)%kT / Kelvin < 10._dp ) then
-          call die('transiesta electronic temperature *must* &
+          call die('TranSiesta electronic temperature *must* &
                &be larger than 10 kT')
        end if
        
@@ -1023,7 +1023,7 @@
        end select
     end if
     if ( .not. Calc_Forces ) then
-       write(*,f11) '*** TranSIESTA will NOT update forces ***'
+       write(*,f11) '*** TranSiesta will NOT update forces ***'
     end if
 
     if ( TS_RHOCORR_METHOD == 0 ) then
@@ -1112,17 +1112,17 @@
 
     write(*,'(3a)') repeat('*',24),' Begin: TS CHECKS AND WARNINGS ',repeat('*',24)
     if ( FixSpin .and. (TS_HS_save .or. TSmode) ) then
-       write(*,'(a)') 'Fixed spin not possible with TranSIESTA!'
+       write(*,'(a)') 'Fixed spin not possible with TranSiesta!'
        write(*,'(a)') 'Disable TS.HS.Save or FixSpin'
-       write(*,'(a)') 'Electrodes with fixed spin is not possible with Transiesta!'
-       call die('Fixing spin is not possible in transiesta')
+       write(*,'(a)') 'Electrodes with fixed spin is not possible with TranSiesta!'
+       call die('Fixing spin is not possible in TranSiesta')
     end if
 
     if ( .not. TSmode ) then
        if ( TS_HA == TS_HA_ELEC ) then
-          call die('Hartree potiental cannot use electrodes without transiesta')
+          call die('Hartree potiental cannot use electrodes without TranSiesta')
        else if ( TS_HA == TS_HA_ELEC_BOX ) then
-          call die('Hartree potiental cannot use electrodes without transiesta')
+          call die('Hartree potiental cannot use electrodes without TranSiesta')
        end if
     end if
 
@@ -1148,7 +1148,7 @@
     end if
 
     if ( TS_HA == TS_HA_NONE ) then
-       write(*,'(a)') 'Hartree potiental fix REQUIRED when running transiesta'
+       write(*,'(a)') 'Hartree potiental fix REQUIRED when running TranSiesta'
        err = .true.
     end if
 
@@ -1310,17 +1310,17 @@
     call contour_nEq_warnings()
     
     if ( .not. Calc_Forces ) then
-       write(*,f11) '***       TranSIESTA will NOT update forces       ***'
+       write(*,f11) '***       TranSiesta will NOT update forces       ***'
        write(*,f11) '*** ALL FORCES AFTER TRANSIESTA HAS RUN ARE WRONG ***'
        if ( Nmove > 0 ) then
-          write(*,'(a)')'Relaxation with transiesta *REQUIRES* an update of &
+          write(*,'(a)')'Relaxation with TranSiesta *REQUIRES* an update of &
                &the energy density matrix. Will continue at your request.'
           err = .true.
        end if
     end if
 
     if ( Nmove > 0 .and. .not. all(Elecs(:)%DM_update > 0) ) then
-       write(*,'(a)') 'transiesta relaxation is only allowed if you also &
+       write(*,'(a)') 'TranSiesta relaxation is only allowed if you also &
             &update, at least, the cross terms, please set: &
             &TS.Elecs.DM.Update [cross-terms|all]'
        err = .true.
@@ -1382,7 +1382,7 @@
 
     if ( ts_tidx < 1 ) then
 
-       write(*,'(a)') '*** TranSIESTA semi-infinite directions are individual ***'
+       write(*,'(a)') '*** TranSiesta semi-infinite directions are individual ***'
        write(*,'(a)') '*** It is heavily adviced to have any electrodes with no &
             &periodicity'
        write(*,'(a)') '    in the transverse directions be located as far from any &
@@ -1626,7 +1626,7 @@
        write(*,'(tr19,a)') 'TRANSIESTA REPORTED ERRORS'
        write(*,'(tr18,a)') repeat('*',30)
 
-       call die('One or more errors have occured doing transiesta &
+       call die('One or more errors have occured doing TranSiesta &
             &initialization, check the output')
     end if
     
@@ -1654,13 +1654,13 @@
 
     if ( onlyS .or. .not. TSmode ) return
 
-    write(*,'(/,a,/)') '>>> Transiesta block information for FDF-file START <<<'
+    write(*,'(/,a,/)') '>>> TranSiesta block information for FDF-file START <<<'
     
     call print_mus_block( 'TS' , N_mu , mus)
     
     call print_contour_block( 'TS' , IsVolt )
     
-    write(*,'(/,a,/)') '>>> Transiesta block information for FDF-file END <<<'
+    write(*,'(/,a,/)') '>>> TranSiesta block information for FDF-file END <<<'
 
     ! write out the contour
     call io_contour(IsVolt, mus, slabel)

=== modified file 'Src/m_ts_pivot.F90'
--- Src/m_ts_pivot.F90	2018-03-26 19:52:10 +0000
+++ Src/m_ts_pivot.F90	2018-05-04 19:08:45 +0000
@@ -314,18 +314,41 @@
        call sp_pvt(n,tmp_Sp,r_pvt, PVT_PCG, c_pvt, start = r_Els, only_sub = .not. lextend )
 
 #ifdef SIESTA__METIS
-    else if ( str_contain(pvt_str,'metis+priority') ) then
-       str_tmp = trim(str_tmp)//'+metis+priority'
+    else if ( str_contain(pvt_str,'metis+priority') .or. &
+        str_contain(pvt_str,'nodend+priority') ) then
+       str_tmp = trim(str_tmp)//'+nodend+priority'
 
-       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS, c_pvt, &
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_NODEND, c_pvt, &
             priority = priority%r)
-    else if ( str_contain(pvt_str,'metis') ) then
-       str_tmp = trim(str_tmp)//'+metis'
-
-       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS, c_pvt)
+    else if ( str_contain(pvt_str,'metis') .or. &
+        str_contain(pvt_str,'nodend') ) then
+       str_tmp = trim(str_tmp)//'+nodend'
+
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_NODEND, c_pvt)
+
+     else if ( str_contain(pvt_str,'partgraphkway+priority') ) then
+       str_tmp = trim(str_tmp)//'+partgraphkway+priority'
+
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_PARTGRAPHKWAY, c_pvt, &
+           priority = priority%r)
+       
+    else if ( str_contain(pvt_str,'partgraphkway') ) then
+       str_tmp = trim(str_tmp)//'+partgraphkway'
+
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_PARTGRAPHKWAY, c_pvt)
+
+    else if ( str_contain(pvt_str,'partgraphrecursive+priority') ) then
+       str_tmp = trim(str_tmp)//'+partgraphrecursive+priority'
+
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_PARTGRAPHRECURSIVE, c_pvt, &
+           priority = priority%r)
+       
+    else if ( str_contain(pvt_str,'partgraphrecursive') ) then
+       str_tmp = trim(str_tmp)//'+partgraphrecursive'
+
+       call sp_pvt(n,tmp_Sp,r_pvt, PVT_METIS_PARTGRAPHRECURSIVE, c_pvt)
 
 #endif
-
     else if ( str_contain(pvt_str,'CG') .or. pvt_option == PVT_NONE ) then
 
        is_rev = str_contain(pvt_str, 'rev-')

=== modified file 'Src/m_ts_tri_init.F90'
--- Src/m_ts_tri_init.F90	2018-04-07 18:57:43 +0000
+++ Src/m_ts_tri_init.F90	2018-05-04 19:08:45 +0000
@@ -631,6 +631,69 @@
        call tri(r_El)
     end if
 
+#ifdef SIESTA__METIS
+    fmethod = trim(corb)//'+NodeND+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'NodeND+priority'
+    call sp_pvt(n,tmpSp2,r_tmp, PVT_METIS_NODEND, sub = full, priority = priority%r)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+    
+    fmethod = trim(corb)//'+rev-NodeND+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'rev-NodeND+priority'
+    call rgn_reverse(r_tmp)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+    
+    fmethod = trim(corb)//'+PartGraphKway+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'PartGraphKway+priority'
+    call sp_pvt(n,tmpSp2,r_tmp, PVT_METIS_PARTGRAPHKWAY, sub = full, priority = priority%r)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+    
+    fmethod = trim(corb)//'+rev-PartGraphKway+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'rev-PartGraphKway+priority'
+    call rgn_reverse(r_tmp)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+
+    fmethod = trim(corb)//'+PartGraphRecursive+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'PartGraphRecursive+priority'
+    call sp_pvt(n,tmpSp2,r_tmp, PVT_METIS_PARTGRAPHRECURSIVE, sub = full, priority = priority%r)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+    
+    fmethod = trim(corb)//'+rev-PartGraphRecursive+priority'
+    if ( IONode ) write(*,fmt) trim(corb),'rev-PartGraphRecursive+priority'
+    call rgn_reverse(r_tmp)
+    if ( orb_atom == 1 ) then
+      call tri(r_tmp)
+    else
+      call rgn_atom2orb(r_tmp,na_u,lasto,r_El)
+      call tri(r_El)
+    end if
+
+#endif
+    
     end do orb_atom_switch
 
     call rgn_delete(r_tmp,r_El,full,priority)

=== modified file 'Util/SpPivot/pvtsp.F90'
--- Util/SpPivot/pvtsp.F90	2018-04-06 17:22:45 +0000
+++ Util/SpPivot/pvtsp.F90	2018-05-04 19:08:45 +0000
@@ -144,9 +144,15 @@
         fmethod = 'Scramble'
         
 #ifdef SIESTA__METIS
-      case ( '-metis' )
-        method = PVT_METIS
-        fmethod = 'Metis'
+      case ( '-metis', '-nodend' )
+        method = PVT_METIS_NODEND
+        fmethod = 'metis-NodeND'
+      case ( '-partgraphkway' )
+        method = PVT_METIS_PARTGRAPHKWAY
+        fmethod = 'metis-PartGraphKway'
+      case ( '-partgraphrecursive' )
+        method = PVT_METIS_PARTGRAPHRECURSIVE
+        fmethod = 'metis-PartGraphRecursive'
 #endif
 
       case ( '-w', '-weight' )
@@ -388,9 +394,9 @@
 
   subroutine help()
     character(len=20), parameter :: gf = '(tr3,a,'':'',/,tr8,a)'
-    character(len=10), parameter :: nf = '(tr8,a)'
+    character(len=*), parameter :: nf = '(tr8,a)'
 
-    character(len=10), parameter :: fm = '(tr11,a)'
+    character(len=*), parameter :: fm = '(tr11,a18,": ",a)'
 
     write(*,'(a)') 'The following options are available for pvtsp:'
     write(*,'(a)') 
@@ -401,18 +407,20 @@
     write(*,gf) '--metis-stdout','make METIS output (on STDOUT)'
     write(*,gf) '--pvt <method>','pivot according to a specific method'
     write(*,nf) '--<method> can be one of the following:'
-    write(*,fm) '      cm: Cuthill-Mckee'
-    write(*,fm) '  rev-cm: reverse Cuthill-Mckee'
-    write(*,fm) '     gps: Gibbs-Poole-Stockmeyer'
-    write(*,fm) ' rev-gps: reverse Gibbs-Poole-Stockmeyer'
-    write(*,fm) '     pcg: Peripheral connectivity graph'
-    write(*,fm) ' rev-pcg: reverse Peripheral connectivity graph'
-    write(*,fm) '    ggps: General Gibbs-Poole-Stockmeyer'
-    write(*,fm) 'rev-ggps: reverse General Gibbs-Poole-Stockmeyer'
+    write(*,fm) 'cm', 'Cuthill-Mckee'
+    write(*,fm) 'rev-cm', 'reverse Cuthill-Mckee'
+    write(*,fm) 'gps', 'Gibbs-Poole-Stockmeyer'
+    write(*,fm) 'rev-gps', 'reverse Gibbs-Poole-Stockmeyer'
+    write(*,fm) 'pcg', 'Peripheral connectivity graph'
+    write(*,fm) 'rev-pcg', 'reverse Peripheral connectivity graph'
+    write(*,fm) 'ggps', 'General Gibbs-Poole-Stockmeyer'
+    write(*,fm) 'rev-ggps', 'reverse General Gibbs-Poole-Stockmeyer'
 #ifdef SIESTA__METIS
-    write(*,fm) '   metis: Metis pivoting'
+    write(*,fm) 'nodend', 'Metis NodeND pivoting'
+    write(*,fm) 'partgraphkway', 'Metis PartGraphKway pivoting'
+    write(*,fm) 'partgraphrecursive', 'Metis PartGraphRecursive pivoting'
 #endif
-    write(*,fm) 'scramble: Scramble the sparsity pattern'
+    write(*,fm) 'scramble', 'Scramble the sparsity pattern'
     write(*,'(a)')
     write(*,gf) '--unit-cell|-uc','convert to unit-cell sparsity pattern'
     write(*,gf) '--a|-a <i>','use only ith supercell as connectivity graph (in A direction)'

=== modified file 'Util/TS/TBtrans/Makefile'
--- Util/TS/TBtrans/Makefile	2018-04-07 19:13:42 +0000
+++ Util/TS/TBtrans/Makefile	2018-05-04 19:08:45 +0000
@@ -445,9 +445,8 @@
 denmatlomem.o: alloc.o globalise.o onmod.o precision.o
 densematrix.o: alloc.o precision.o
 detover.o: alloc.o linpack.o parallel.o parallelsubs.o precision.o
-dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o
-dfscf.o: m_spin.o mesh.o meshdscf.o meshphi.o parallel.o parallelsubs.o
-dfscf.o: precision.o
+dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o mesh.o
+dfscf.o: meshdscf.o meshphi.o parallel.o parallelsubs.o precision.o
 dhscf.o: alloc.o atmfuncs.o bsc_xcmod.o cellxc_mod.o delk.o dfscf.o
 dhscf.o: doping_uniform.o files.o forhar.o iogrid_netcdf.o local_sys.o
 dhscf.o: m_charge_add.o m_efield.o m_hartree_add.o m_iorho.o m_mesh_node.o
@@ -655,11 +654,11 @@
 m_new_dm.o: class_Fstack_Pair_Geometry_SpData2D.o class_Geometry.o
 m_new_dm.o: class_OrbitalDistribution.o class_Pair_Geometry_SpData2D.o
 m_new_dm.o: class_SpData2D.o class_Sparsity.o files.o local_sys.o m_energies.o
-m_new_dm.o: m_handle_sparse.o m_iodm.o m_mixing.o m_mixing_scf.o m_mpi_utils.o
-m_new_dm.o: m_spin.o m_spin.o m_steps.o m_svd.o m_ts_electype.o
-m_new_dm.o: m_ts_global_vars.o m_ts_iodm.o m_ts_method.o m_ts_options.o
-m_new_dm.o: parallel.o parsing.o precision.o restructSpData2D.o siesta_geom.o
-m_new_dm.o: siesta_options.o sparse_matrices.o units.o
+m_new_dm.o: m_handle_sparse.o m_iodm.o m_mixing.o m_mixing_scf.o m_spin.o
+m_new_dm.o: m_spin.o m_steps.o m_svd.o m_ts_electype.o m_ts_global_vars.o
+m_new_dm.o: m_ts_iodm.o m_ts_method.o m_ts_options.o parallel.o parsing.o
+m_new_dm.o: precision.o restructSpData2D.o siesta_geom.o siesta_options.o
+m_new_dm.o: sparse_matrices.o units.o
 m_noccbands.o: alloc.o atmfuncs.o local_sys.o m_spin.o parallel.o precision.o
 m_noccbands.o: siesta_geom.o
 m_options.o: precision.o
@@ -733,13 +732,14 @@
 m_ts_contour_eq.o: m_gauss_fermi_17.o m_gauss_fermi_18.o m_gauss_fermi_19.o
 m_ts_contour_eq.o: m_gauss_fermi_20.o m_gauss_fermi_22.o m_gauss_fermi_24.o
 m_ts_contour_eq.o: m_gauss_fermi_26.o m_gauss_fermi_28.o m_gauss_fermi_30.o
-m_ts_contour_eq.o: m_gauss_fermi_inf.o m_gauss_quad.o m_integrate.o m_ts_aux.o
-m_ts_contour_eq.o: m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
+m_ts_contour_eq.o: m_gauss_fermi_inf.o m_gauss_quad.o m_integrate.o m_io.o
+m_ts_contour_eq.o: m_ts_aux.o m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
 m_ts_contour_eq.o: m_ts_io_contour.o m_ts_io_ctype.o parallel.o precision.o
 m_ts_contour_eq.o: units.o
-m_ts_contour_neq.o: m_gauss_quad.o m_integrate.o m_ts_aux.o m_ts_cctype.o
-m_ts_contour_neq.o: m_ts_chem_pot.o m_ts_electype.o m_ts_io_contour.o
-m_ts_contour_neq.o: m_ts_io_ctype.o parallel.o precision.o units.o
+m_ts_contour_neq.o: m_gauss_quad.o m_integrate.o m_io.o m_ts_aux.o
+m_ts_contour_neq.o: m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
+m_ts_contour_neq.o: m_ts_io_contour.o m_ts_io_ctype.o parallel.o precision.o
+m_ts_contour_neq.o: units.o
 m_ts_debug.o: class_Sparsity.o class_TriMat.o geom_helper.o parallel.o
 m_ts_debug.o: precision.o
 m_ts_dm_update.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
@@ -782,7 +782,7 @@
 m_ts_io.o: class_SpData2D.o class_Sparsity.o geom_helper.o local_sys.o m_io_s.o
 m_ts_io.o: m_ncdf_io.o m_os.o m_sparse.o parallel.o precision.o
 m_ts_io_contour.o: precision.o units.o
-m_ts_io_ctype.o: parallel.o precision.o units.o
+m_ts_io_ctype.o: m_io.o parallel.o precision.o units.o
 m_ts_iodm.o: class_OrbitalDistribution.o class_SpData2D.o class_Sparsity.o
 m_ts_iodm.o: m_io_s.o m_os.o parallel.o precision.o
 m_ts_kpoints.o: files.o find_kgrid.o kpoint_grid.o local_sys.o
@@ -904,7 +904,7 @@
 mulliken.o: siesta_cml.o
 naefs.o: atmfuncs.o mneighb.o new_matel.o precision.o
 new_matel.o: alloc.o errorf.o interpolation.o local_sys.o matel_registry.o
-new_matel.o: parallel.o precision.o radfft.o spher_harm.o
+new_matel.o: parallel.o precision.o radfft.o siesta_options.o spher_harm.o
 nlefsm.o: alloc.o atm_types.o atmfuncs.o atomlist.o chemical.o mneighb.o
 nlefsm.o: new_matel.o parallel.o parallelsubs.o precision.o
 normalize_dm.o: atomlist.o local_sys.o m_mpi_utils.o m_spin.o parallel.o
@@ -990,9 +990,10 @@
 save_density_matrix.o: m_ts_global_vars.o m_ts_iodm.o m_ts_options.o
 save_density_matrix.o: precision.o siesta_options.o sparse_matrices.o
 savepsi.o: alloc.o parallel.o parallelsubs.o precision.o
-scfconvergence_test.o: ldau.o ldau_specs.o m_convergence.o m_energies.o
-scfconvergence_test.o: m_wallclock.o parallel.o precision.o siesta_cml.o
-scfconvergence_test.o: siesta_options.o units.o write_subs.o
+scfconvergence_test.o: atomlist.o ldau.o ldau_specs.o m_convergence.o
+scfconvergence_test.o: m_energies.o m_spin.o m_wallclock.o parallel.o
+scfconvergence_test.o: precision.o siesta_cml.o siesta_geom.o siesta_options.o
+scfconvergence_test.o: sparse_matrices.o units.o write_subs.o
 schecomm.o: alloc.o
 setatomnodes.o: alloc.o local_sys.o parallel.o precision.o spatial.o
 setspatial.o: alloc.o parallel.o precision.o spatial.o
@@ -1064,10 +1065,10 @@
 siesta_move.o: atomlist.o broyden_optim.o cell_broyden_optim.o
 siesta_move.o: cell_fire_optim.o dynamics.o fire_optim.o flook_siesta.o ioxv.o
 siesta_move.o: local_sys.o m_check_walltime.o m_energies.o m_exp_coord.o
-siesta_move.o: m_forces.o m_kinetic.o m_steps.o m_stress.o m_target_stress.o
-siesta_move.o: parallel.o siesta_cml.o siesta_dicts.o siesta_geom.o
-siesta_move.o: siesta_master.o siesta_options.o units.o write_subs.o
-siesta_move.o: zm_broyden_optim.o zm_fire_optim.o zmatrix.o
+siesta_move.o: m_forces.o m_kinetic.o m_mpi_utils.o m_steps.o m_stress.o
+siesta_move.o: m_target_stress.o parallel.o siesta_cml.o siesta_dicts.o
+siesta_move.o: siesta_geom.o siesta_master.o siesta_options.o units.o
+siesta_move.o: write_subs.o zm_broyden_optim.o zm_fire_optim.o zmatrix.o
 sparse_matrices.o: alloc.o class_Fstack_Pair_Geometry_SpData2D.o
 sparse_matrices.o: class_OrbitalDistribution.o class_SpData1D.o
 sparse_matrices.o: class_SpData2D.o class_Sparsity.o precision.o
@@ -1140,7 +1141,7 @@
 zmatrix.o: alloc.o local_sys.o m_cell.o parallel.o precision.o siesta_geom.o
 zmatrix.o: units.o
 local_sys.o: parallel.o
-m_tbt_contour.o: m_gauss_quad.o m_integrate.o m_tbt_save.o m_ts_aux.o
+m_tbt_contour.o: m_gauss_quad.o m_integrate.o m_io.o m_tbt_save.o m_ts_aux.o
 m_tbt_contour.o: m_ts_cctype.o m_ts_chem_pot.o m_ts_electype.o
 m_tbt_contour.o: m_ts_io_contour.o m_ts_io_ctype.o parallel.o precision.o
 m_tbt_contour.o: units.o
@@ -1177,15 +1178,16 @@
 m_tbt_regions.o: class_OrbitalDistribution.o class_Sparsity.o
 m_tbt_regions.o: create_Sparsity_SC.o create_Sparsity_Union.o fdf_extra.o
 m_tbt_regions.o: files.o geom_helper.o intrinsic_missing.o m_char.o m_pivot.o
-m_tbt_regions.o: m_pivot_methods.o m_region.o m_sparsity_handling.o
-m_tbt_regions.o: m_ts_debug.o m_ts_electype.o m_ts_method.o m_ts_pivot.o
-m_tbt_regions.o: m_ts_sparse.o m_verbosity.o parallel.o precision.o
+m_tbt_regions.o: m_pivot_methods.o m_region.o m_sparsity_handling.o m_tbt_dH.o
+m_tbt_regions.o: m_tbt_delta.o m_ts_debug.o m_ts_electype.o m_ts_method.o
+m_tbt_regions.o: m_ts_pivot.o m_ts_sparse.o m_verbosity.o parallel.o
+m_tbt_regions.o: precision.o
 m_tbt_save.o: class_OrbitalDistribution.o class_SpData1D.o class_Sparsity.o
 m_tbt_save.o: files.o m_interpolate.o m_ncdf_io.o m_os.o m_region.o m_tbt_hs.o
 m_tbt_save.o: m_ts_electype.o m_verbosity.o parallel.o precision.o timestamp.o
 m_tbt_save.o: units.o
 m_tbt_sigma_save.o: m_os.o m_region.o m_tbt_hs.o m_tbt_save.o m_ts_electype.o
-m_tbt_sigma_save.o: parallel.o timestamp.o units.o
+m_tbt_sigma_save.o: parallel.o precision.o timestamp.o units.o
 m_tbt_sparse_helper.o: class_OrbitalDistribution.o class_SpData1D.o
 m_tbt_sparse_helper.o: class_Sparsity.o geom_helper.o intrinsic_missing.o
 m_tbt_sparse_helper.o: m_region.o m_tbt_kregions.o m_ts_electype.o

=== modified file 'Util/TS/TBtrans/m_tbt_contour.F90'
--- Util/TS/TBtrans/m_tbt_contour.F90	2018-03-22 10:43:42 +0000
+++ Util/TS/TBtrans/m_tbt_contour.F90	2018-05-04 19:08:45 +0000
@@ -239,15 +239,19 @@
     type(ts_cw), intent(inout) :: c
     real(dp), intent(in) :: Eta
 
-    if ( leqi(c%c_io%part,'line') ) then
-       
-       call contour_line(c,Eta)
-       
+    if ( leqi(c%c_io%part,'user') ) then
+       
+      call contour_file(c,Eta)
+      
+    else if ( leqi(c%c_io%part,'line') ) then
+      
+      call contour_line(c,Eta)
+      
     else
-       
-       call neq_die('Unrecognized contour type for &
-            &tbtrans, MUST be a line part.')
-       
+      
+      call neq_die('Unrecognized contour type for &
+          &tbtrans, MUST be a line part.')
+      
     end if
     
   end subroutine setup_tbt_contour
@@ -308,7 +312,14 @@
        end if
 
        call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp)
-       
+
+    case ( CC_USER )
+
+      call contour_file(c,Eta)
+
+      deallocate(ce, cw)
+      return
+      
     case default
 
        call die('Could not determine the line-integral')
@@ -322,6 +333,108 @@
     
   end subroutine contour_line
 
+  ! This routine will read the contour points from a file
+  subroutine contour_file(c,Eta)
+    use m_io, only: io_assign, io_close
+    use fdf, only: fdf_convfac
+
+    type(ts_cw), intent(inout) :: c
+    ! The lifting into the complex plane
+    real(dp), intent(in) :: Eta
+
+    integer :: iu, iostat, ne
+    logical :: exist
+    complex(dp) :: E , W
+    real(dp) :: rE, iE, rW, iW, conv
+    character(len=512) :: file, line
+    character(len=16) :: unit
+    
+    ! The contour type contains the file name in:
+    !  c%c_io%cN (weirdly enough)
+    file = c%c_io%cN
+
+    call io_assign(iu)
+    
+    ! Open the contour file
+    inquire(file=trim(file), exist=exist)
+    if ( .not. exist ) then
+      call die('The file: '//trim(file)//' could not be found to read contour points!')
+    end if
+    
+    ! Open the file
+    open(iu, file=trim(file), form='formatted', status='old')
+    
+    ne = 0
+    ! The default unit is eV.
+    ! On every line an optional unit-specificer may be used to specify the
+    ! subsequent lines units (until another unit is specified)
+    conv = fdf_convfac('eV', 'Ry')
+    do
+      ! Now we have the line
+      read(iu, '(a)', iostat=iostat) line
+      if ( iostat /= 0 ) exit
+      if ( len_trim(line) == 0 ) cycle
+      line = trim(adjustl(line))
+      if ( line(1:1) == '#' ) cycle
+
+      ! We have a line with energy and weight
+      ne = ne + 1
+      ! There are three optional ways of reading this
+      ! 1.  ReE ImE, ReW ImW [unit]
+      read(line, *, iostat=iostat) rE, iE, rW, iW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW, iW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 2.  ReE ImE, ReW [unit]
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, iE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, iE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE,iE) * conv
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      ! 3.  ReE , ReW [unit]
+      iE = Eta
+      iW = 0._dp
+      read(line, *, iostat=iostat) rE, rW, unit
+      if ( iostat == 0 ) then
+        conv = fdf_convfac(unit, 'Ry')
+      else
+        read(line, *, iostat=iostat) rE, rW
+      end if
+      if ( iostat == 0 ) then
+        c%c(ne) = dcmplx(rE * conv,iE)
+        c%w(ne,1) = dcmplx(rW,iW) * conv
+        cycle
+      end if
+
+      call die('Contour file: '//trim(file)//' is not formatted correctly. &
+          &Please read the documentation!')
+
+    end do
+
+    call io_close(iu)
+
+    if ( c%c_io%N /= ne ) then
+      call die('Error in reading the contour points from file: '//trim(file))
+    end if
+    
+  end subroutine contour_file
+
   function TBT_E(id,step) result(c)
     integer, intent(in) :: id
     integer, intent(in), optional :: step

=== modified file 'Util/TS/TBtrans/m_tbt_options.F90'
--- Util/TS/TBtrans/m_tbt_options.F90	2018-03-22 10:43:42 +0000
+++ Util/TS/TBtrans/m_tbt_options.F90	2018-05-04 19:08:45 +0000
@@ -602,6 +602,26 @@
             &apply.','Set TBT.DOS.A T to calculate orbital currents.'
     end if
 
+    ! Options for density-matrix calculations
+    ltmp = fdf_get('TBT.DM.Gf', .false.)
+    if ( ltmp .and. ('DOS-Gf'.in.save_DATA)) then
+      save_DATA = save_DATA // ('DM-Gf'.kv.1)
+    else if ( ltmp .and. IONode ) then
+      write(*,'(2(/,a))')'WARNING: Will not calculate the density matrix (Gf), &
+          &the Green function DOS needs to be calculated for this to &
+          &apply.','Set TBT.DOS.Gf T to calculate density matrix (Gf).'
+    end if
+
+    ltmp = fdf_get('TBT.DM.A', .false.)
+    if ( ltmp .and. ('DOS-A'.in.save_DATA)) then
+       save_DATA = save_DATA // ('DM-A'.kv.1)
+    else if ( ltmp .and. IONode ) then
+       write(*,'(2(/,a))')'WARNING: Will not calculate the density matrix (A), &
+            &the spectral function DOS needs to be calculated for this to &
+            &apply.','Set TBT.DOS.A T to calculate density matrix (A).'
+    end if
+
+    
     ! Options for COOP and COHP curves.
     ! These are orbital (energy) populations that can be used to determine the
     ! bonding nature of the material.
@@ -619,7 +639,7 @@
        save_DATA = save_DATA // ('COOP-A'.kv.1)
     else if ( ltmp .and. IONode ) then
        write(*,'(2(/,a))')'WARNING: Will not calculate the COOP (A) curve, &
-            &the Green function DOS needs to be calculated for this to &
+            &the spectral function DOS needs to be calculated for this to &
             &apply.','Set TBT.DOS.A T to calculate COOP (A) curves.'
     end if
 
@@ -637,7 +657,7 @@
        save_DATA = save_DATA // ('COHP-A'.kv.1)
     else if ( ltmp .and. IONode ) then
        write(*,'(2(/,a))')'WARNING: Will not calculate the COHP (A) curve, &
-            &the Green function DOS needs to be calculated for this to &
+            &the spectral function DOS needs to be calculated for this to &
             &apply.','Set TBT.DOS.A T to calculate COHP (A) curves.'
     end if
 
@@ -749,6 +769,10 @@
     end if
     write(*,f1) 'Saving bond currents (orb-orb)',('orb-current'.in.save_DATA)
 
+    ! DM
+    write(*,f1) 'Saving DM from Green function',('DM-Gf'.in.save_DATA)
+    write(*,f1) 'Saving DM from spectral functions',('DM-A'.in.save_DATA)
+
     ! COOP/COHP curves
     write(*,f1) 'Saving COOP from Green function',('COOP-Gf'.in.save_DATA)
     write(*,f1) 'Saving COOP from spectral functions',('COOP-A'.in.save_DATA)
@@ -874,7 +898,7 @@
        has = has .or. ('COOP-A' .in. save_DATA)
        has = has .or. ('COHP-Gf' .in. save_DATA)
        has = has .or. ('COHP-A' .in. save_DATA)
-       if ( IONode .and. ltmp .and. has) then
+       if ( IONode .and. ltmp .and. has ) then
           write(*,'(a,/,a)') 'WARNING: k-averaging COOP/COHP with &
                &time-reversal symmetry will not reproduce','the correct &
                &populations. Set TBT.Symmetry.TimeReversal F'

=== modified file 'Util/TS/TBtrans/m_tbt_proj.F90'
--- Util/TS/TBtrans/m_tbt_proj.F90	2018-03-22 10:43:42 +0000
+++ Util/TS/TBtrans/m_tbt_proj.F90	2018-05-04 19:08:45 +0000
@@ -1306,7 +1306,7 @@
 
   ! Initialize the TBT.Proj.nc file
   subroutine init_proj_save( fname, TSHS , r, btd, ispin, N_Elec, Elecs, &
-       nkpt, kpt, wkpt, NE , a_Dev, a_Buf, sp_dev_sc, save_DATA )
+       nkpt, kpt, wkpt, NE , Eta, a_Dev, a_Buf, sp_dev_sc, save_DATA )
 
     use parallel, only : Node, Nodes, IONode
     use units, only: eV
@@ -1341,6 +1341,7 @@
     integer, intent(in) :: N_Elec
     type(Elec), intent(in) :: Elecs(N_Elec)
     integer, intent(in) :: nkpt, NE
+    real(dp), intent(in) :: Eta
     real(dp), intent(in) :: kpt(3,nkpt), wkpt(nkpt)
     type(tRgn), intent(in) :: a_Dev
     type(tRgn), intent(in) :: a_Buf
@@ -1684,8 +1685,11 @@
 #else
     dic = dic//('info'.kv.'Energy')//('unit'.kv.'Ry')
 #endif
-    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), &
-         atts = dic)
+    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), atts = dic)
+
+    dic = dic//('info'.kv.'Imaginary part for device')
+    call ncdf_def_var(ncdf,'eta',NF90_DOUBLE,(/'one'/), atts = dic)
+
     call delete(dic)
 
     call ncdf_put_var(ncdf,'nsc',TSHS%nsc)
@@ -1716,6 +1720,8 @@
     call ncdf_put_var(ncdf,'wkpt',wkpt)
     deallocate(rv)
 
+    call ncdf_put_var(ncdf,'eta',Eta)
+
     if ( 'proj-orb-current' .in. save_DATA ) then
 
        ! In case we need to save the device sparsity pattern

=== modified file 'Util/TS/TBtrans/m_tbt_regions.F90'
--- Util/TS/TBtrans/m_tbt_regions.F90	2018-04-07 18:57:43 +0000
+++ Util/TS/TBtrans/m_tbt_regions.F90	2018-05-04 19:08:45 +0000
@@ -773,7 +773,11 @@
 #ifdef MPI
     use mpi_siesta, only : MPI_Comm_Self
 #endif
-    use m_sparsity_handling, only : Sp_retain_region, Sp_sort
+    use m_sparsity_handling, only : Sp_retain_region, Sp_sort, Sp_union
+#ifdef NCDF_4
+    use m_tbt_delta, only: read_delta_Sp
+    use m_tbt_dH, only: use_dH, dH
+#endif
 
     type(Sparsity), intent(inout) :: sp
     type(dict), intent(in) :: save_DATA
@@ -781,6 +785,7 @@
     type(OrbitalDistribution) :: fdit
 
     integer :: no_u
+    type(Sparsity) :: sp_dH
 #endif
 
     ! Make sure to initialize the device region
@@ -789,6 +794,7 @@
 #ifdef NCDF_4
     if ( ('orb-current' .in. save_DATA) .or. &
          ('proj-orb-current' .in. save_DATA) .or. &
+         ('DM-Gf' .in. save_DATA) .or. ('DM-A' .in. save_DATA) .or. &
          ('COOP-Gf' .in. save_DATA) .or. ('COHP-Gf' .in. save_DATA) .or. &
          ('COOP-A' .in. save_DATA) .or. ('COHP-A' .in. save_DATA) ) then
 
@@ -798,8 +804,15 @@
 #else
        call newDistribution(no_u,-1           ,fdit,name='TBT-fake dist')
 #endif
-
        call Sp_retain_region(fdit,sp,r_oDev,sp_dev_sc)
+       ! Note that the delta-Sigma is not necessary because
+       ! the self-energy does not add to bond-currents, etc.
+       if ( use_dH ) then
+         call read_delta_Sp(dH,no_u,sp_dH)
+         call Sp_retain_region(fdit,sp_dH,r_oDev,sp_dH)
+         call Sp_union(fdit,sp_dev_sc,sp_dH,sp_dev_sc)
+         call delete(sp_dH)
+       end if
        call Sp_sort(sp_dev_sc)
        call delete(fdit)
 

=== modified file 'Util/TS/TBtrans/m_tbt_save.F90'
--- Util/TS/TBtrans/m_tbt_save.F90	2018-04-09 12:39:44 +0000
+++ Util/TS/TBtrans/m_tbt_save.F90	2018-05-04 19:08:45 +0000
@@ -311,7 +311,7 @@
 
   subroutine init_cdf_save(fname,TSHS,r,btd,ispin, &
       N_Elec, Elecs, rEl, btd_El, &
-      nkpt, kpt, wkpt, NE, &
+      nkpt, kpt, wkpt, NE, Eta, &
       a_Dev, a_Buf, sp_dev_sc, &
       save_DATA )
 
@@ -349,6 +349,7 @@
     integer, intent(in) :: nkpt
     real(dp), intent(in), target :: kpt(3,nkpt), wkpt(nkpt)
     integer, intent(in) :: NE
+    real(dp), intent(in) :: Eta
     type(tRgn), intent(in) :: a_Dev
     ! In case the system has some buffer atoms.
     type(tRgn), intent(in) :: a_Buf
@@ -362,7 +363,7 @@
     type(dict) :: dic
     logical :: exist, sme, isGamma
     integer :: iEl, jEl, i, nnzs_dev, N_eigen
-    integer :: prec_DOS, prec_T, prec_Teig, prec_J, prec_COOP
+    integer :: prec_DOS, prec_T, prec_Teig, prec_J, prec_COOP, prec_DM
     type(OrbitalDistribution) :: fdit
     real(dp) :: mem
     character(len=2) :: unit
@@ -383,6 +384,7 @@
     call tbt_cdf_precision('T.Eig','single',prec_Teig)
     call tbt_cdf_precision('Current','single',prec_J)
     call tbt_cdf_precision('COOP','single',prec_COOP)
+    call tbt_cdf_precision('DM','single',prec_DM)
 
     isGamma = all(TSHS%nsc(:) == 1)
 
@@ -579,6 +581,7 @@
     dic = dic // ('info'.kv.'Blocks in BTD for the pivot table')
     call ncdf_def_var(ncdf,'btd',NF90_INT,(/'n_btd'/), &
          atts = dic)
+    mem = mem + calc_mem(NF90_INT, btd%n)
 
     dic = dic // ('info'.kv.'Index of device atoms')
     call ncdf_def_var(ncdf,'a_dev',NF90_INT,(/'na_d'/), &
@@ -613,11 +616,15 @@
 #else
     dic = dic//('info'.kv.'Energy')//('unit'.kv.'Ry')
 #endif
-    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), &
-         atts = dic, chunks = (/1/) )
-    call delete(dic)
+    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), atts = dic)
     mem = mem + calc_mem(NF90_DOUBLE, NE)
 
+    dic = dic//('info'.kv.'Imaginary part for device')
+    call ncdf_def_var(ncdf,'eta',NF90_DOUBLE,(/'one'/), atts = dic)
+
+    ! Clean-up dictionary
+    call delete(dic)
+
     call ncdf_put_var(ncdf,'nsc',TSHS%nsc)
     call ncdf_put_var(ncdf,'isc_off',TSHS%isc_off)
     call ncdf_put_var(ncdf,'pivot',r%r)
@@ -646,11 +653,15 @@
     call ncdf_put_var(ncdf,'wkpt',wkpt)
     deallocate(r2)
 
+    call ncdf_put_var(ncdf,'eta',Eta)
+
     sme = 'orb-current' .in. save_DATA
     sme = sme .or. ('COOP-Gf' .in. save_DATA)
     sme = sme .or. ('COOP-A' .in. save_DATA)
     sme = sme .or. ('COHP-Gf' .in. save_DATA)
     sme = sme .or. ('COHP-A' .in. save_DATA)
+    sme = sme .or. ('DM-Gf' .in. save_DATA)
+    sme = sme .or. ('DM-A' .in. save_DATA)
     if ( sme ) then
        
        ! In case we need to save the device sparsity pattern
@@ -683,6 +694,12 @@
 
     end if
 
+    if ( 'DM-Gf' .in. save_DATA ) then
+       dic = dic // ('info'.kv.'Green function density matrix')
+       call ncdf_def_var(ncdf,'DM',prec_DM,(/'nnzs','ne  ','nkpt'/), &
+           atts = dic , chunks = (/nnzs_dev/) , compress_lvl=cmp_lvl)
+       mem = mem + calc_mem(prec_DM, nnzs_dev, NE, nkpt)
+    end if
     if ( 'COOP-Gf' .in. save_DATA ) then
        dic = dic // ('info'.kv.'Crystal orbital overlap population')//('unit'.kv.'1/Ry')
        call ncdf_def_var(ncdf,'COOP',prec_COOP,(/'nnzs','ne  ','nkpt'/), &
@@ -747,15 +764,18 @@
           dic = ('info'.kv.'Last orbitals of the equivalent atom')
           call ncdf_def_var(grp,'lasto',NF90_INT,(/'na_u'/), &
                atts = dic)
+          mem = mem + calc_mem(NF90_INT, Elecs(iEl)%na_u)
+          
           dic = dic//('info'.kv.'Bulk transmission')
           call ncdf_def_var(grp,'T',prec_T,(/'ne  ','nkpt'/), &
-              atts = dic)
+               atts = dic)
           mem = mem + calc_mem(prec_T, NE, nkpt)
 
           dic = dic//('info'.kv.'Unit cell')
           dic = dic//('unit'.kv.'Bohr')
           call ncdf_def_var(grp,'cell',NF90_DOUBLE,(/'xyz','xyz'/), &
                atts = dic)
+          
           dic = dic//('info'.kv.'Atomic coordinates')
           call ncdf_def_var(grp,'xa',NF90_DOUBLE,(/'xyz ','na_u'/), &
                atts = dic , chunks = (/3, Elecs(iEl)%na_u/) )
@@ -773,7 +793,7 @@
           mem = mem + calc_mem(prec_DOS, Elecs(iEl)%no_u, NE, nkpt)
 
        end if
-
+       call delete(dic)
        
        ! Now we will only add information that is calculated
        if ( iEl == N_Elec ) then
@@ -782,6 +802,12 @@
                ('T-all'.nin. save_DATA) ) cycle
        end if
 
+       if ( 'DM-A' .in. save_DATA ) then
+          dic = dic // ('info'.kv.'Spectral function density matrix')
+          call ncdf_def_var(grp,'DM',prec_DM,(/'nnzs','ne  ','nkpt'/), &
+              atts = dic , chunks = (/nnzs_dev/) , compress_lvl=cmp_lvl)
+          mem = mem + calc_mem(prec_DM, nnzs_dev, NE, nkpt)
+       end if
 
        if ( 'DOS-A' .in. save_DATA ) then
           dic = dic//('info'.kv.'Spectral function density of states')// &
@@ -1162,7 +1188,7 @@
 #endif
 
 #ifdef TBTRANS_TIMING
-    call timer('cdf-w-DTJ',1)
+    call timer('cdf-w-DOS-T',1)
 #endif
 
     NDOS = size(DOS,dim=1)
@@ -1308,7 +1334,7 @@
 #endif
 
 #ifdef TBTRANS_TIMING
-    call timer('cdf-w-DTJ',2)
+    call timer('cdf-w-DOS-T',2)
 #endif
 
   end subroutine state_cdf_save

=== modified file 'Util/TS/TBtrans/m_tbt_sigma_save.F90'
--- Util/TS/TBtrans/m_tbt_sigma_save.F90	2018-04-09 12:39:44 +0000
+++ Util/TS/TBtrans/m_tbt_sigma_save.F90	2018-05-04 19:08:45 +0000
@@ -139,7 +139,7 @@
 
   ! Save the self-energies of the electrodes and
   subroutine init_Sigma_save(fname, TSHS, r, btd, ispin, N_Elec, Elecs, &
-       nkpt, kpt, wkpt, NE, &
+       nkpt, kpt, wkpt, NE, Eta, &
        a_Dev, a_Buf)
 
     use parallel, only : IONode
@@ -168,6 +168,8 @@
     integer, intent(in) :: nkpt
     real(dp), intent(in) :: kpt(3,nkpt), wkpt(nkpt)
     integer, intent(in) :: NE
+    real(dp), intent(in) :: Eta
+
     ! Device atoms
     type(tRgn), intent(in) :: a_Dev
     ! Buffer atoms
@@ -379,8 +381,11 @@
 #else
     dic = dic//('info'.kv.'Energy')//('unit'.kv.'Ry')
 #endif
-    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), &
-         atts = dic)
+    call ncdf_def_var(ncdf,'E',NF90_DOUBLE,(/'ne'/), atts = dic)
+
+    dic = dic//('info'.kv.'Imaginary part for device')
+    call ncdf_def_var(ncdf,'eta',NF90_DOUBLE,(/'one'/), atts = dic)
+
     call delete(dic)
 
     call ncdf_put_var(ncdf,'pivot',r%r)
@@ -412,6 +417,8 @@
     deallocate(r2)
     mem = mem + calc_mem(NF90_DOUBLE, 4, nkpt)
 
+    call ncdf_put_var(ncdf,'eta',Eta)
+
     do iEl = 1 , N_Elec
 
        call ncdf_def_grp(ncdf,trim(Elecs(iEl)%name),grp)

=== modified file 'Util/TS/TBtrans/m_tbt_tri_scat.F90'
--- Util/TS/TBtrans/m_tbt_tri_scat.F90	2017-12-07 21:13:32 +0000
+++ Util/TS/TBtrans/m_tbt_tri_scat.F90	2018-05-04 19:08:45 +0000
@@ -47,6 +47,7 @@
   public :: GF_COHP_add_dH, A_COHP_add_dH
   public :: orb_current
   public :: orb_current_add_dH
+  public :: GF_DM, A_DM
 #endif
 
   ! Used for BLAS calls (local variables)
@@ -68,18 +69,19 @@
   ! This routine utilizes the sparse matrix as a loop, instead of looping
   ! all BTD matrix elements.
   ! This turns out to be much faster for (at least tight-binding calculations).
-  subroutine GF_DOS(r,Gf_tri,work_tri,S_1D,pvt,DOS)
+  subroutine GF_DOS(r,Gfd_tri,Gfo_tri,S_1D,pvt,DOS)
     use class_Sparsity
     use class_zSpData1D
     
     type(tRgn), intent(in) :: r
-    type(zTriMat), intent(inout) :: Gf_tri, work_tri
+    type(zTriMat), intent(inout) :: Gfd_tri, Gfo_tri
     type(zSpData1D), intent(inout) :: S_1D ! (transposed S(k))
     type(tRgn), intent(in) :: pvt
     real(dp), intent(out) :: DOS(r%n)
 
     type(Sparsity), pointer :: sp
     complex(dp), pointer :: S(:)
+    complex(dp), pointer :: Gfd(:), Gfo(:)
     complex(dp) :: GfGfd
     integer, pointer :: ncol(:), l_ptr(:), l_col(:)
     integer :: np, n, no_o, no_i
@@ -87,29 +89,29 @@
     real(dp) :: lDOS
 
 #ifdef TBTRANS_TIMING
-    call timer('GF-DOS',1)
+    call timer('Gf-DOS',1)
 #endif
 
-    np = parts(Gf_tri)
+    np = parts(Gfd_tri)
     
     ! First calculate all off-diagonal green-function elements
-    no_o = nrows_g(Gf_tri,1)
-    no_i = nrows_g(Gf_tri,2)
+    no_o = nrows_g(Gfd_tri,1)
+    no_i = nrows_g(Gfd_tri,2)
     call calc(2,1)
     do n = 2, np - 1
-       no_o = nrows_g(Gf_tri,n)
-       no_i = nrows_g(Gf_tri,n + 1)
+       no_o = nrows_g(Gfd_tri,n)
+       no_i = nrows_g(Gfd_tri,n + 1)
        call calc(n+1,n)
-       no_i = nrows_g(Gf_tri,n - 1)
+       no_i = nrows_g(Gfd_tri,n - 1)
        call calc(n-1,n)
     end do
-    no_o = nrows_g(Gf_tri,np)
-    no_i = nrows_g(Gf_tri,np-1)
+    no_o = nrows_g(Gfd_tri,np)
+    no_i = nrows_g(Gfd_tri,np-1)
     call calc(np-1,np)
 
     ! At this point we have calculated all Green function matrices
-    ! All diagonal elements are in Gf_tri,
-    ! all off-diagonal elements are in work_tri
+    ! All diagonal elements are in Gfd_tri,
+    ! all off-diagonal elements are in Gfo_tri
 
     ! The DOS per orbital is calculated like this (.=matrix multiplication):
     !   DOS(io) = - Im[ (Gf-Gf^\dagger) . S ](io,io) / Pi
@@ -126,6 +128,9 @@
     S => val(S_1D)
     call attach(sp,n_col=ncol,list_ptr=l_ptr,list_col=l_col)
 
+    Gfd => val(Gfd_tri)
+    Gfo => val(Gfo_tri)
+
 !$OMP parallel do default(shared), private(br,io,lDOS,ind,bc,GfGfd)
     do br = 1, r%n
        io = r%r(br)
@@ -135,7 +140,7 @@
        do ind = l_ptr(io) + 1, l_ptr(io) + ncol(io)
           bc = pvt%r(l_col(ind))
           if ( bc > 0 ) then
-             call index_GfGfd(br, bc, GfGfd)
+             call calc_GfGfd(br, bc, GfGfd)
              lDOS = lDOS + dimag( GfGfd * S(ind) )
           end if
        end do
@@ -146,7 +151,7 @@
 !$OMP end parallel do
 
 #ifdef TBTRANS_TIMING
-    call timer('GF-DOS',2)
+    call timer('Gf-DOS',2)
 #endif
 
   contains
@@ -155,9 +160,9 @@
       integer, intent(in) :: m,n
       complex(dp), pointer :: Gf(:), Mnn(:), XY(:)
 
-      XY  => val(Gf_tri,m,n)
-      Mnn => val(Gf_tri,n,n)
-      Gf  => val(work_tri,m,n)
+      XY  => val(Gfd_tri,m,n)
+      Mnn => val(Gfd_tri,n,n)
+      Gf  => val(Gfo_tri,m,n)
       
       ! We need to calculate the 
       ! Mnm1n/Mnp1n Green's function
@@ -171,29 +176,26 @@
       
     end subroutine calc
 
-    subroutine index_GfGfd(br, bc, G)
+    subroutine calc_GfGfd(br, bc, G)
       integer, intent(in) :: br, bc
-      complex(dp), intent(out) :: G
-      complex(dp), pointer :: Gf(:)
-      integer :: p_r, i_r, p_c, i_c
-
-      call part_index(work_tri, br, p_r, i_r)
-      call part_index(work_tri, bc, p_c, i_c)
-
+      complex(dp), intent(inout) :: G
+      integer :: p_r, i_r, p_c, i_c, i
+
+      call part_index(Gfo_tri, br, p_r, i_r)
+      call part_index(Gfo_tri, bc, p_c, i_c)
+      
       if ( p_r == p_c ) then
-         Gf => val(Gf_tri, p_r, p_c)
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfd(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        G = G - conjg(Gfd(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
       else
-         Gf => val(work_tri, p_r, p_c)
-      end if
-      G = Gf(i_r + (i_c-1) * work_tri%data%tri_nrows(p_r))
-
-      ! Immediate subtract G^\dagger element
-      if ( p_r /= p_c ) then
-         Gf => val(work_tri, p_c, p_r)
-      end if
-      G = G - conjg(Gf(i_c + (i_r-1) * work_tri%data%tri_nrows(p_c)))
-
-    end subroutine index_GfGfd
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfo(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        i = index_block(Gfo_tri, p_c, p_r)
+        G = G - conjg(Gfo(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
+      end if
+
+    end subroutine calc_GfGfd
 
   end subroutine GF_DOS
 
@@ -267,7 +269,7 @@
   ! This routine utilizes the sparse matrix as a loop, instead of looping
   ! all BTD matrix elements.
   ! This turns out to be much faster for (at least tight-binding calculations).
-  subroutine GF_COP(r,Gfd_tri,Gfo_tri,pvt,sp,M,sc_off,k,COP)
+  subroutine GF_COP(r,Gfd_tri,Gfo_tri,pvt,sp,M,sc_off,k,ph,COP)
     use class_Sparsity
     use class_dSpData1D
     use intrinsic_missing, only : SFIND
@@ -280,18 +282,19 @@
     real(dp), intent(in) :: M(:) ! S for COOP, H for COHP
     real(dp), intent(in) :: sc_off(:,:)
     real(dp), intent(in) :: k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(dSpData1D), intent(inout) :: COP ! COOP or COHP
 
     type(Sparsity), pointer :: c_sp
     real(dp), pointer :: C(:)
+    complex(dp), pointer :: Gfd(:), Gfo(:)
     complex(dp) :: GfGfd
-    complex(dp), allocatable :: ph(:)
     integer, pointer :: ncol(:), l_ptr(:), l_col(:)
     integer, pointer :: cncol(:), cptr(:), ccol(:), c_col(:)
     integer :: no_u, br, io, ind, iind, bc
 
 #ifdef TBTRANS_TIMING
-    call timer('GF-COP',1)
+    call timer('Gf-COP',1)
 #endif
 
 #ifdef TBT_PHONON
@@ -317,7 +320,6 @@
     ! Since we have to do Gf.S we simply
     ! create S(-k) (which is S^T)
     ! and thus get the correct values.
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do io = 1 , size(sc_off, dim=2)
        ph(io-1) = cdexp(dcmplx(0._dp, - &
             k(1) * sc_off(1,io) - &
@@ -331,6 +333,9 @@
     C => val(COP)
     call attach(c_sp, n_col=cncol, list_ptr=cptr, list_col=ccol)
 
+    Gfd => val(Gfd_tri)
+    Gfo => val(Gfo_tri)
+
 !$OMP parallel default(shared), private(br,io,c_col,ind,iind,bc,GfGfd)
 
 !$OMP workshare
@@ -357,7 +362,7 @@
 
           ! COOP(iind) = - Im[ (G(io,jo) - G^\dagger(io,jo)) * S(jo,io) ] / 2Pi
           bc = pvt%r(ucorb(l_col(ind),no_u)) ! pivoted orbital index in tri-diagonal matrix
-          call index_GfGfd(br, bc, GfGfd)
+          call calc_GfGfd(br, bc, GfGfd)
           
           C(iind) = -aimag( GfGfd * M(ind) * ph( (l_col(ind)-1)/no_u ))
 
@@ -367,42 +372,36 @@
 !$OMP end do
 !$OMP end parallel
 
-    ! Clean-up phases
-    deallocate(ph)
-    
 #ifdef TBTRANS_TIMING
-    call timer('GF-COP',2)
+    call timer('Gf-COP',2)
 #endif
 
   contains
 
-    subroutine index_GfGfd(br, bc, G)
+    subroutine calc_GfGfd(br, bc, G)
       integer, intent(in) :: br, bc
-      complex(dp), intent(out) :: G
-      complex(dp), pointer :: Gf(:)
-      integer :: p_r, i_r, p_c, i_c
-
+      complex(dp), intent(inout) :: G
+      integer :: p_r, i_r, p_c, i_c, i
+      
       call part_index(Gfo_tri, br, p_r, i_r)
       call part_index(Gfo_tri, bc, p_c, i_c)
-
+      
       if ( p_r == p_c ) then
-         Gf => val(Gfd_tri, p_r, p_c)
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfd(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        G = G - conjg(Gfd(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
       else
-         Gf => val(Gfo_tri, p_r, p_c)
-      end if
-      G = Gf(i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
-
-      ! Immediate subtract G^\dagger element
-      if ( p_r /= p_c ) then
-         Gf => val(Gfo_tri, p_c, p_r)
-      end if
-      G = G - conjg(Gf(i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
-
-    end subroutine index_GfGfd
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfo(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        i = index_block(Gfo_tri, p_c, p_r)
+        G = G - conjg(Gfo(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
+      end if
+
+    end subroutine calc_GfGfd
 
   end subroutine GF_COP
 
-  subroutine Gf_COHP_add_dH(dH_1D,sc_off,k,Gfd_tri,Gfo_tri,r,COHP,pvt)
+  subroutine Gf_COHP_add_dH(dH_1D,sc_off,k,ph,Gfd_tri,Gfo_tri,r,COHP,pvt)
 
     use class_Sparsity
     use class_zSpData1D
@@ -412,6 +411,7 @@
 
     type(zSpData1D), intent(in) :: dH_1D
     real(dp), intent(in) :: sc_off(:,:), k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(zTriMat), intent(inout) :: Gfd_tri, Gfo_tri
     type(tRgn), intent(in) :: r
     type(dSpData1D), intent(inout) :: COHP
@@ -424,7 +424,7 @@
     integer, pointer :: cncol(:), cptr(:), ccol(:)
     integer, pointer :: l_ncol(:), l_ptr(:), l_col(:), col(:)
 
-    complex(dp), allocatable :: ph(:)
+    complex(dp), pointer :: Gfd(:), Gfo(:)
     complex(dp) :: GfGfd
     real(dp), pointer :: C(:)
     integer :: no_u, br, io, jo, i, ind, iind
@@ -444,14 +444,16 @@
     call attach(c_sp, n_col=cncol, list_ptr=cptr, list_col=ccol)
 
     ! Create the phases
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do i = 1 , size(sc_off, dim=2)
-       ph(i-1) = cdexp(dcmplx(0._dp, &
-            k(1) * sc_off(1,i) + &
-            k(2) * sc_off(2,i) + &
+       ph(i-1) = cdexp(dcmplx(0._dp, - &
+            k(1) * sc_off(1,i) - &
+            k(2) * sc_off(2,i) - &
             k(3) * sc_off(3,i))) / (2._dp * Pi)
     end do
 
+    Gfd => val(Gfd_tri)
+    Gfo => val(Gfo_tri)
+
 !$OMP parallel do default(shared), private(br,io,iind,jo,ind,col,GfGfd)
     do br = 1, r%n
        io = r%r(br)
@@ -476,7 +478,7 @@
           
              if ( ind > l_ptr(jo) ) then
                 
-                call index_GfGfd(br, pvt%r(jo), GfGfd)
+                call calc_GfGfd(br, pvt%r(jo), GfGfd)
                 ! COHP(iind) += - Im[ (G(io,jo) - G^\dagger(io,jo)) * dH(jo,io)] / 2Pi
                 C(iind) = C(iind) &
                      - aimag( GfGfd * dH(ind) * ph( (l_col(ind)-1)/no_u ))
@@ -490,37 +492,32 @@
     end do
 !$OMP end parallel do
 
-    deallocate(ph)
-
 #ifdef TBTRANS_TIMING
     call timer('COHP-Gf-dH',2)
 #endif
 
   contains
 
-    subroutine index_GfGfd(br, bc, G)
+    subroutine calc_GfGfd(br, bc, G)
       integer, intent(in) :: br, bc
-      complex(dp), intent(out) :: G
-      complex(dp), pointer :: Gf(:)
-      integer :: p_r, i_r, p_c, i_c
+      complex(dp), intent(inout) :: G
+      integer :: p_r, i_r, p_c, i_c, i
 
       call part_index(Gfo_tri, br, p_r, i_r)
       call part_index(Gfo_tri, bc, p_c, i_c)
-
+      
       if ( p_r == p_c ) then
-         Gf => val(Gfd_tri, p_r, p_c)
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfd(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        G = G - conjg(Gfd(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
       else
-         Gf => val(Gfo_tri, p_r, p_c)
-      end if
-      G = Gf(i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
-
-      ! Immediate subtract G^\dagger element
-      if ( p_r /= p_c ) then
-         Gf => val(Gfo_tri, p_c, p_r)
-      end if
-      G = G - conjg(Gf(i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
-
-    end subroutine index_GfGfd
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfo(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        i = index_block(Gfo_tri, p_c, p_r)
+        G = G - conjg(Gfo(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
+      end if
+
+    end subroutine calc_GfGfd
 
     function TO(io) result(jo)
       integer, intent(in) :: io
@@ -553,7 +550,7 @@
   ! This routine utilizes the sparse matrix as a loop, instead of looping
   ! all BTD matrix elements.
   ! This turns out to be much faster for (at least tight-binding calculations).
-  subroutine A_COP(r,A_tri,pvt,sp,M,sc_off,k,COP)
+  subroutine A_COP(r,A_tri,pvt,sp,M,sc_off,k,ph,COP)
     use class_Sparsity
     use class_dSpData1D
     use intrinsic_missing, only : SFIND
@@ -566,12 +563,12 @@
     real(dp), intent(in) :: M(:) ! S for COOP, H for COHP
     real(dp), intent(in) :: sc_off(:,:)
     real(dp), intent(in) :: k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(dSpData1D), intent(inout) :: COP ! COOP or COHP
 
     type(Sparsity), pointer :: c_sp
     real(dp), pointer :: C(:)
     complex(dp), pointer :: A(:)
-    complex(dp), allocatable :: ph(:)
     integer, pointer :: ncol(:), l_ptr(:), l_col(:)
     integer, pointer :: cncol(:), cptr(:), ccol(:), c_col(:)
     integer :: no_u, br, io, ind, iind, bc
@@ -590,7 +587,7 @@
 
     ! The COOP calculation can be written as
     !
-    !   COOP(io,jo) = - Im{ A(io,jo) * S(jo,io) * e^(ik.R) } / 2Pi
+    !   COOP(io,jo) = Re{ A(io,jo) * S(jo,io) * e^(ik.R) } / 2Pi
     ! Here we want:
     !   ADOS(io) = \sum_jo COOP(io,jo)
     ! since we know that COOP(io,jo) is the io -> jo ADOS.
@@ -601,7 +598,6 @@
     ! Since we have to do A.S we simply
     ! create the S(-k) (which is S^T)
     ! and thus get the correct values.
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do io = 1 , size(sc_off, dim=2)
        ph(io-1) = cdexp(dcmplx(0._dp, - &
             k(1) * sc_off(1,io) - &
@@ -653,16 +649,13 @@
 !$OMP end do
 !$OMP end parallel
 
-    ! Clean-up phases
-    deallocate(ph)
-    
 #ifdef TBTRANS_TIMING
     call timer('A-COP',2)
 #endif
 
   end subroutine A_COP
 
-  subroutine A_COHP_add_dH(dH_1D,sc_off,k,A_tri,r,COHP,pvt)
+  subroutine A_COHP_add_dH(dH_1D,sc_off,k,ph,A_tri,r,COHP,pvt)
 
     use class_Sparsity
     use class_zSpData1D
@@ -672,6 +665,7 @@
 
     type(zSpData1D), intent(in) :: dH_1D
     real(dp), intent(in) :: sc_off(:,:), k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(zTriMat), intent(inout) :: A_tri
     type(tRgn), intent(in) :: r
     type(dSpData1D), intent(inout) :: COHP
@@ -684,7 +678,6 @@
     integer, pointer :: cncol(:), cptr(:), ccol(:)
     integer, pointer :: l_ncol(:), l_ptr(:), l_col(:), col(:)
 
-    complex(dp), allocatable :: ph(:)
     complex(dp), pointer :: A(:)
     real(dp), pointer :: C(:)
     integer :: no_u, iu, io, i, ind, iind, jo, iA
@@ -704,9 +697,8 @@
     call attach(c_sp, n_col=cncol, list_ptr=cptr, list_col=ccol)
 
     ! Create the phases
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do i = 1 , size(sc_off, dim=2)
-       ph(i-1) = cdexp(dcmplx(0._dp, &
+       ph(i-1) = cdexp(dcmplx(0._dp, + &
             k(1) * sc_off(1,i) + &
             k(2) * sc_off(2,i) + &
             k(3) * sc_off(3,i))) / (2._dp * Pi)
@@ -750,8 +742,6 @@
     end do
 !$OMP end parallel do
 
-    deallocate(ph)
-
 #ifdef TBTRANS_TIMING
     call timer('COHP-A-dH',2)
 #endif
@@ -1568,7 +1558,7 @@
 
 
 #ifdef NCDF_4
-  subroutine orb_current(sp,H,S,sc_off,k,cE,A_tri,r,orb_J,pvt)
+  subroutine orb_current(sp,H,S,sc_off,k,ph,cE,A_tri,r,orb_J,pvt)
 
     use class_Sparsity
     use class_zSpData1D
@@ -1582,6 +1572,7 @@
     ! We require that the input Hamiltonian is Hermitian
     real(dp), intent(in) :: H(:), S(:), sc_off(:,:)
     real(dp), intent(in) :: k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(ts_c_idx) :: cE
     type(zTriMat), intent(inout) :: A_tri
     ! The region that specifies the size of orb_J
@@ -1594,7 +1585,6 @@
     integer, pointer :: i_ncol(:), i_ptr(:), i_col(:), icol(:)
     integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
 
-    complex(dp), allocatable :: ph(:)
     complex(dp), pointer :: A(:)
     complex(dp) :: Hi
     real(dp), pointer :: J(:)
@@ -1615,12 +1605,9 @@
     J    => val (orb_J)
     call attach(i_sp, n_col=i_ncol, list_ptr=i_ptr, list_col=i_col)
 
-    ! We do not initialize J as every entry is overwritten
-
     ! Create the phases
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do i = 1 , size(sc_off, dim=2)
-       ph(i-1) = cdexp(dcmplx(0._dp, &
+       ph(i-1) = cdexp(dcmplx(0._dp, + &
             k(1) * sc_off(1,i) + &
             k(2) * sc_off(2,i) + &
             k(3) * sc_off(3,i)))
@@ -1685,15 +1672,13 @@
 !$OMP end do
 !$OMP end parallel
 
-    deallocate(ph)
-
 #ifdef TBTRANS_TIMING
     call timer('orb-current',2)
 #endif
 
   end subroutine orb_current
   
-  subroutine orb_current_add_dH(dH_1D,sc_off,k,A_tri,r,orb_J,pvt)
+  subroutine orb_current_add_dH(dH_1D,sc_off,k,ph,A_tri,r,orb_J,pvt)
 
     use class_Sparsity
     use class_zSpData1D
@@ -1703,6 +1688,7 @@
 
     type(zSpData1D), intent(in) :: dH_1D
     real(dp), intent(in) :: sc_off(:,:), k(3)
+    complex(dp), intent(inout) :: ph(0:)
     type(zTriMat), intent(inout) :: A_tri
     ! The region that specifies the size of orb_J
     type(tRgn), intent(in) :: r
@@ -1716,7 +1702,6 @@
     integer, pointer :: i_ncol(:), i_ptr(:), i_col(:)
     integer, pointer :: l_ncol(:), l_ptr(:), l_col(:), col(:)
 
-    complex(dp), allocatable :: ph(:)
     complex(dp) :: p
     complex(dp), pointer :: A(:)
     real(dp), pointer :: J(:)
@@ -1737,9 +1722,8 @@
     call attach(i_sp, n_col=i_ncol, list_ptr=i_ptr, list_col=i_col)
 
     ! Create the phases
-    allocate( ph(0:size(sc_off,dim=2)-1) )
     do i = 1 , size(sc_off, dim=2)
-       ph(i-1) = cdexp(dcmplx(0._dp, &
+       ph(i-1) = cdexp(dcmplx(0._dp, + &
             k(1) * sc_off(1,i) + &
             k(2) * sc_off(2,i) + &
             k(3) * sc_off(3,i)))
@@ -1817,8 +1801,6 @@
     end do
 !$OMP end parallel do
 
-    deallocate(ph)
-
 #ifdef TBTRANS_TIMING
     call timer('orb-current-dH',2)
 #endif
@@ -1851,6 +1833,183 @@
     end function TO
     
   end subroutine orb_current_add_dH
+
+
+  subroutine GF_DM(sc_off,k,ph,Gfd_tri,Gfo_tri,r,pvt,spDM)
+
+    use class_Sparsity
+    use class_dSpData1D
+    use geom_helper,       only : UCORB
+
+    real(dp), intent(in) :: sc_off(:,:)
+    real(dp), intent(in) :: k(3)
+    complex(dp), intent(inout) :: ph(0:)
+    type(zTriMat), intent(inout) :: Gfd_tri, Gfo_tri
+    ! The region that specifies the size of spDM
+    type(tRgn), intent(in) :: r
+    ! The pivoting region that transfers r%r(iu) to io
+    type(tRgn), intent(in) :: pvt
+    type(dSpData1D), intent(inout) :: spDM
+
+    integer, pointer :: ncol(:), l_ptr(:), l_col(:)
+
+    type(Sparsity), pointer :: sp
+    real(dp), pointer :: DM(:)
+    complex(dp), pointer :: Gfd(:), Gfo(:)
+    complex(dp) :: GfGfd
+
+    integer :: no_u, iu, io, ind, ju
+
+#ifdef TBTRANS_TIMING
+    call timer('Gf-DM',1)
+#endif
+
+    sp => spar(spDM)
+    DM => val(spDM)
+    call attach(sp, nrows_g=no_u, n_col=ncol, list_ptr=l_ptr, list_col=l_col)
+
+    ! Create the phases
+    do io = 1 , size(sc_off, dim=2)
+       ph(io-1) = cdexp(dcmplx(0._dp, - &
+            k(1) * sc_off(1,io) - &
+            k(2) * sc_off(2,io) - &
+            k(3) * sc_off(3,io))) / (2._dp * Pi)
+    end do
+
+    Gfd => val(Gfd_tri)
+    Gfo => val(Gfo_tri)
+    
+!$OMP parallel default(shared), private(iu,io,ind,ju,GfGfd)
+
+    ! we need this in case the device region gets enlarged due to dH
+!$OMP workshare
+    DM(:) = 0._dp
+!$OMP end workshare
+    
+!$OMP do
+    do iu = 1, r%n
+       io = r%r(iu)
+
+#ifndef TS_NOCHECKS
+       if ( ncol(io) == 0 ) call die('Gf_DM: DM has zero columns &
+            &for at least one row')
+#endif
+
+       ! Loop on DM entries here...
+       do ind = l_ptr(io) + 1 , l_ptr(io) + ncol(io)
+
+         ju = pvt%r(ucorb(l_col(ind), no_u))
+         call calc_GfGfd(iu, ju, GfGfd)
+         DM(ind) = - aimag( GfGfd * ph((l_col(ind) - 1) / no_u) )
+
+       end do
+    end do
+!$OMP end do
+!$OMP end parallel
+
+#ifdef TBTRANS_TIMING
+    call timer('Gf-DM',2)
+#endif
+    
+  contains
+    
+    subroutine calc_GfGfd(br, bc, G)
+      integer, intent(in) :: br, bc
+      complex(dp), intent(inout) :: G
+      integer :: p_r, i_r, p_c, i_c, i
+
+      call part_index(Gfo_tri, br, p_r, i_r)
+      call part_index(Gfo_tri, bc, p_c, i_c)
+      
+      if ( p_r == p_c ) then
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfd(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        G = G - conjg(Gfd(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
+      else
+        i = index_block(Gfo_tri, p_r, p_c)
+        G = Gfo(i + i_r + (i_c-1) * Gfo_tri%data%tri_nrows(p_r))
+        i = index_block(Gfo_tri, p_c, p_r)
+        G = G - conjg(Gfo(i + i_c + (i_r-1) * Gfo_tri%data%tri_nrows(p_c)))
+      end if
+
+    end subroutine calc_GfGfd
+    
+  end subroutine GF_DM
+
+  subroutine A_DM(sc_off,k,ph,A_tri,r,pvt,spDM)
+
+    use class_Sparsity
+    use class_dSpData1D
+    use geom_helper,       only : UCORB
+
+    real(dp), intent(in) :: sc_off(:,:)
+    real(dp), intent(in) :: k(3)
+    complex(dp), intent(inout) :: ph(0:)
+    type(zTriMat), intent(inout) :: A_tri
+    ! The region that specifies the size of spDM
+    type(tRgn), intent(in) :: r
+    ! The pivoting region that transfers r%r(iu) to io
+    type(tRgn), intent(in) :: pvt
+    type(dSpData1D), intent(inout) :: spDM
+
+    integer, pointer :: ncol(:), l_ptr(:), l_col(:)
+
+    type(Sparsity), pointer :: sp
+    complex(dp), pointer :: A(:)
+    real(dp), pointer :: DM(:)
+    integer :: no_u, iu, io, ind, ju
+
+#ifdef TBTRANS_TIMING
+    call timer('A-DM',1)
+#endif
+
+    sp => spar(spDM)
+    DM => val(spDM)
+    call attach(sp, nrows_g=no_u, n_col=ncol, list_ptr=l_ptr, list_col=l_col)
+
+    ! Create the phases
+    do io = 1 , size(sc_off, dim=2)
+       ph(io-1) = cdexp(dcmplx(0._dp, - &
+            k(1) * sc_off(1,io) - &
+            k(2) * sc_off(2,io) - &
+            k(3) * sc_off(3,io))) / (2._dp * Pi)
+    end do
+
+    A => val(A_tri)
+!$OMP parallel default(shared), private(iu,io,ind,ju)
+
+    ! we need this in case the device region gets enlarged due to dH
+!$OMP workshare
+    DM(:) = 0._dp
+!$OMP end workshare
+    
+!$OMP do
+    do iu = 1, r%n
+       io = r%r(iu)
+
+#ifndef TS_NOCHECKS
+       if ( ncol(io) == 0 ) call die('A_DM: DM has zero columns &
+            &for at least one row')
+#endif
+
+       ! Loop on DM entries here...
+       do ind = l_ptr(io) + 1 , l_ptr(io) + ncol(io)
+
+         ju = pvt%r(ucorb(l_col(ind), no_u))
+         ju = index(A_tri, iu, ju)
+         DM(ind) = real(A(ju) * ph((l_col(ind) - 1) / no_u), dp)
+
+       end do
+    end do
+!$OMP end do
+!$OMP end parallel
+
+#ifdef TBTRANS_TIMING
+    call timer('A-DM',2)
+#endif
+
+  end subroutine A_DM
+
 #endif
 
   subroutine insert_Self_energy(n1,n2,M,r,El,off1,off2)
@@ -1968,202 +2127,4 @@
 
   end subroutine insert_Self_energy_Dev
 
-
-#ifdef THESE_HAVE_BEEN_SUPERSEEDED
-  
-  ! A simple routine to calculate the DOS
-  ! from a partially calculated GF
-  ! When entering this routine Gf_tri
-  ! should contain:
-  ! all GF_nn
-  ! all Yn/Bn-1 and all Xn/Cn+1
-  ! This lets us calculate all entries
-  subroutine GF_DOS_loop_BTD(r,Gf_tri,S_1D,DOS,nwork,work)
-    use intrinsic_missing, only : SFIND
-    use class_Sparsity
-    use class_zSpData1D
-
-    type(tRgn), intent(in) :: r
-    type(zTriMat), intent(inout) :: Gf_tri
-    type(zSpData1D), intent(inout) :: S_1D
-    real(dp), intent(out) :: DOS(r%n)
-    integer, intent(in) :: nwork
-    complex(dp), intent(inout), target :: work(nwork)
-
-    type(Sparsity), pointer :: sp
-    complex(dp), pointer :: S(:), Gf(:), Mnn(:), XY(:)
-    integer, pointer :: ncol(:), l_ptr(:), l_col(:), lcol(:)
-    integer :: off1, off2, n, in
-    integer :: jo, i, j, no_o, no_i, ind, np
-    real(dp) :: lDOS
-
-#ifdef TBTRANS_TIMING
-    call timer('GF-DOS',1)
-#endif
-
-    S  => val(S_1D)
-    sp => spar(S_1D)
-    call attach(sp,n_col=ncol,list_ptr=l_ptr,list_col=l_col)
-    
-    ! Initialize DOS to 0
-!$OMP parallel workshare default(shared)
-    DOS(:) = 0._dp
-!$OMP end parallel workshare
-
-    off2 = 0
-    np = parts(Gf_tri)
-    do n = 1 , np
-
-       no_o = nrows_g(Gf_tri,n)
-
-       do in = max(1,n-1) , min(n+1,np)
-
-          no_i = nrows_g(Gf_tri,in)
-
-          if ( in < n ) then
-             off1 = off2 - no_i
-          else if ( n < in ) then
-             off1 = off2 + no_o
-          else
-             off1 = off2
-          end if
-
-          if ( in == n ) then
-             ! Retrieve the central part of the
-             ! matrix
-             Gf => val(Gf_tri,n,n)
-
-          else
-
-             XY  => val(Gf_tri,in,n)
-             Mnn => val(Gf_tri,n,n)
-
-             Gf  => work(1:no_o*no_i)
-
-             ! We need to calculate the 
-             ! Mnm1n/Mnp1n Green's function
-#ifdef USE_GEMM3M
-             call zgemm3m( &
-#else
-             call zgemm( &
-#endif
-                  'N','N',no_i,no_o,no_o, &
-                  zm1, XY,no_i, Mnn,no_o,z0, Gf,no_i)
-
-          end if
-
-!$OMP parallel do default(shared), private(j,jo,ind,i,lcol,lDOS)
-          do j = 1 , no_o
-             jo = r%r(off2+j)
-             lcol => l_col(l_ptr(jo)+1:l_ptr(jo)+ncol(jo))
-             ! get the equivalent one in the
-             ! overlap matrix
-             ! REMEMBER, S is transposed!
-             ! Hence we do not need conjg :)
-             lDOS = 0._dp
-             do i = 1 , no_i
-                ind = SFIND(lcol,r%r(off1+i))
-                if ( ind == 0 ) cycle
-                ind = l_ptr(jo) + ind
-                lDOS = lDOS - dimag( Gf(j+(i-1)*no_i) * S(ind) )
-             end do
-             DOS(off2+j) = DOS(off2+j) + lDOS / Pi
-          end do
-!$OMP end parallel do
-
-       end do
-
-       ! Update the offset
-       off2 = off2 + no_o
-
-    end do
-
-#ifdef TBTRANS_TIMING
-    call timer('GF-DOS',2)
-#endif
-
-  end subroutine GF_DOS_loop_BTD
-
-  ! A simple routine to calculate the DOS
-  ! from a full calculated spectral function
-  ! This loops over entries in the BTD matrix.
-  subroutine A_DOS_loop_BTD(r,A_tri,S_1D,DOS)
-    use intrinsic_missing, only : SFIND
-    use class_Sparsity
-    use class_zSpData1D
-
-    type(tRgn), intent(in) :: r
-    type(zTriMat), intent(inout) :: A_tri
-    type(zSpData1D), intent(inout) :: S_1D
-    real(dp), intent(out) :: DOS(r%n)
-
-    type(Sparsity), pointer :: sp
-    complex(dp), pointer :: S(:), A(:)
-    integer, pointer :: ncol(:), l_ptr(:), l_col(:), lcol(:)
-    integer :: off1, off2, n, in
-    integer :: jo, i, j, no_o, no_i, ind, np
-    real(dp) :: lDOS
-
-#ifdef TBTRANS_TIMING
-    call timer('A-DOS',1)
-#endif
-
-    S  => val(S_1D)
-    sp => spar(S_1D)
-    call attach(sp,n_col=ncol,list_ptr=l_ptr,list_col=l_col)
-
-    off2 = 0
-    np = parts(A_tri)
-    do n = 1 , np
-
-       no_o = nrows_g(A_tri,n)
-
-       do in = max(1,n-1) , min(n+1,np)
-
-          A => val(A_tri,in,n)
-
-          no_i = nrows_g(A_tri,in)
-
-          if ( in < n ) then
-             off1 = off2 - no_i
-          else if ( n < in ) then
-             off1 = off2 + no_o
-          else
-             off1 = off2
-          end if
-
-!$OMP parallel do default(shared), private(j,jo,ind,i,lcol,lDOS)
-          do j = 1 , no_o
-             jo = r%r(off2+j)
-             lcol => l_col(l_ptr(jo)+1:l_ptr(jo)+ncol(jo))
-             ! get the equivalent one in the
-             ! overlap matrix
-             ! REMEMBER, S is transposed!
-             ! Hence we are doing it correctly
-             lDOS = 0._dp
-             do i = 1 , no_i
-                ind = SFIND(lcol,r%r(off1+i))
-                if ( ind == 0 ) cycle
-                ind = l_ptr(jo) + ind
-                lDOS = lDOS + dreal( A(j+(i-1)*no_i) * S(ind) )
-             end do
-             DOS(off2+j) = DOS(off2+j) + lDOS / (2._dp * Pi)
-          end do
-!$OMP end parallel do
-
-       end do
-
-       ! Update the offset
-       off2 = off2 + no_o
-
-    end do
-
-#ifdef TBTRANS_TIMING
-    call timer('A-DOS',2)
-#endif
-
-  end subroutine A_DOS_loop_BTD
-  
-#endif
-
 end module m_tbt_tri_scat

=== modified file 'Util/TS/TBtrans/m_tbt_trik.F90'
--- Util/TS/TBtrans/m_tbt_trik.F90	2018-03-28 07:50:09 +0000
+++ Util/TS/TBtrans/m_tbt_trik.F90	2018-05-04 19:08:45 +0000
@@ -151,8 +151,12 @@
     real(dp), pointer :: H2D(:,:), S(:), H(:)
     ! To figure out which parts of the tri-diagonal blocks we need
     ! to calculate
-    logical :: calc_DOS_Gf, calc_DOS_A, calc_orb_current
+    logical :: calc_T_Gf, calc_T_out
+    logical :: calc_DOS_Elecs
+    logical :: calc_DOS_Gf, calc_DOS_A
 #ifdef NCDF_4
+    logical :: calc_orb_current
+    logical :: calc_DM_Gf, calc_DM_A
     logical :: calc_COOP_Gf, calc_COOP_A
     logical :: calc_COHP_Gf, calc_COHP_A
 #endif
@@ -188,7 +192,7 @@
     type(tLvlMolEl), pointer :: p_E
     real(dp), allocatable :: pDOS(:,:,:)
     real(dp), allocatable :: bTk(:,:), bTkeig(:,:,:)
-    type(dSpData1D) :: orb_J
+    type(dSpData1D) :: dev_M
 #endif
 ! ************************************************************
 
@@ -196,7 +200,7 @@
     integer :: nGFGGF ! For the triple-product
     complex(dp), pointer :: GFGGF_work(:) => null()
     integer :: ntt_work
-    complex(dp), pointer :: tt_work(:), eig(:)
+    complex(dp), pointer :: tt_work(:), eig(:), phase(:)
 ! ************************************************************
 
 ! ******************* Computational variables ****************
@@ -242,12 +246,17 @@
     real :: last_progress, cur_progress
 ! ************************************************************
 
+    calc_DOS_Elecs = 'DOS-Elecs' .in. save_DATA
+    calc_T_Gf = 'T-Gf' .in. save_DATA
+    calc_T_out = 'T-sum-out' .in. save_DATA
     calc_DOS_Gf = 'DOS-Gf' .in. save_DATA
     calc_DOS_A = 'DOS-A' .in. save_DATA
-    calc_orb_current = 'orb-current' .in. save_DATA
     T_all = 'T-all' .in. save_DATA
     ADOS_all = 'DOS-A-all' .in. save_DATA
 #ifdef NCDF_4
+    calc_orb_current = 'orb-current' .in. save_DATA
+    calc_DM_Gf = 'DM-Gf' .in. save_DATA
+    calc_DM_A = 'DM-A' .in. save_DATA
     calc_COOP_Gf = 'COOP-Gf' .in. save_DATA
     calc_COOP_A = 'COOP-A' .in. save_DATA
     calc_COHP_Gf = 'COHP-Gf' .in. save_DATA
@@ -437,14 +446,6 @@
     end if
 
     
-    ! Note that padding is the extra size to be able to calculate
-    ! the spectral function in the BTD format
-
-    ! We allocate for the matrix
-    call print_memory('LHS',pad_LHS)
-    call newzTriMat(zwork_tri,DevTri%n,DevTri%r,'GFinv', &
-         padding = pad_LHS )
-
     ! The RHS will be the array that retains the 
     ! self-energies and nothing more.
     ! Here we pad with the missing elements to contain
@@ -461,10 +462,25 @@
     ! this is only preferred in tbtrans as there might be 
     ! padding any-way.
     pad_RHS = max(pad_RHS,nGFGGF)
+
+    ! Ensure that at least one of the work arrays has
+    ! elements for the few supercells
+    io = product(TSHS%nsc)
+    if ( pad_LHS < io .and. pad_RHS < io ) then
+      pad_LHS = io
+    end if
     
+    ! Note that padding is the extra size to be able to calculate
+    ! the spectral function in the BTD format
+
+    ! We allocate for the matrix
+    call print_memory('LHS',pad_LHS)
+    call newzTriMat(zwork_tri,DevTri%n,DevTri%r,'GFinv', &
+         padding = pad_LHS )
+
     call print_memory('RHS',pad_RHS)
     call newzTriMat(GF_tri,DevTri%n,DevTri%r,'GF', &
-         padding = pad_RHS )
+        padding = pad_RHS )
 
     ! Create the work-array...
     zwork   => val(zwork_tri,all=.true.)
@@ -472,12 +488,23 @@
     Gfwork  => val(Gf_tri,all=.true.)
     nGfwork =  size(Gfwork)
     if ( nGfwork > nzwork ) then
-       nmaxwork = nGfwork
-       maxwork => Gfwork(:)
-    else
-       nmaxwork = nzwork
-       maxwork => zwork(:)
-    end if
+      nmaxwork = nGfwork
+      maxwork => Gfwork(:)
+    else
+      nmaxwork = nzwork
+      maxwork => zwork(:)
+    end if
+
+    ! Create phase array
+    io = product(TSHS%nsc)
+    if ( pad_LHS >= io ) then
+      phase => zwork(nzwork-io+1:)
+    else if ( pad_RHS >= io ) then
+      phase => Gfwork(nGfwork-io+1:)
+    else
+      call die('Error in coding! Phase size!')
+    end if
+
 
     ! Point the work-array for eigenvalue calculation
     if ( N_eigen > 0 ) then
@@ -601,8 +628,7 @@
        ! The first partition is whether we want
        ! certain quantities from all electrodes
        A_parts(:) = .false.
-       if ( T_all .or. &
-            ('T-sum-out' .in. save_DATA) ) then
+       if ( T_all .or. calc_T_out ) then
 
           ! We need _all_ diagonal blocks of the spectral function
           do iEl = 1 , N_Elec
@@ -651,10 +677,11 @@
 
 #ifdef NCDF_4
     if ( calc_orb_current .or. ('proj-orb-current'.in.save_DATA) .or. &
+         calc_DM_Gf .or. calc_DM_A .or. & 
          calc_COOP_Gf .or. calc_COOP_A .or. & 
          calc_COHP_Gf .or. calc_COHP_A ) then
        
-       call newdSpData1D(sp_dev_sc,fdist,orb_J,name='TBT sparse')
+       call newdSpData1D(sp_dev_sc,fdist,dev_M,name='TBT sparse')
        
     end if
 #endif
@@ -845,7 +872,7 @@
           ! We have reduced the electrode sizes to only one spin-channel
           ! Hence, it will ALWAYS be the first index
           if ( n_k == 0 ) then
-             if ( 'DOS-Elecs' .in. save_DATA ) then
+             if ( calc_DOS_Elecs ) then
                 call read_next_GS(1, ikpt, bkpt, &
                      cE, N_Elec, uGF, Elecs, &
                      nzwork, zwork, .false., forward = .false. , &
@@ -997,7 +1024,7 @@
           if ( .not. only_proj ) then
 #endif
 
-          call timer('DOS-Gf-A-T',1)
+          call timer('analysis',1)
 
           ! We have now calculated all block diagonal entries
           ! of the Green's function.
@@ -1012,19 +1039,23 @@
              end if
              
 #ifdef NCDF_4
+             if ( calc_DM_Gf ) then
+               call Gf_DM(TSHS%sc_off,kpt,phase,Gf_tri,zwork_tri,r_oDev,pvt, dev_M)
+               call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'DM', dev_M)
+             end if
              if ( calc_COOP_Gf ) then
                 call GF_COP(r_oDev,Gf_tri,zwork_tri,pvt, &
-                     TSHS%sp,S,TSHS%sc_off, kpt, orb_J)
-                call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COOP', orb_J)
+                     TSHS%sp,S,TSHS%sc_off, kpt, phase, dev_M)
+                call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COOP', dev_M)
              end if
              if ( calc_COHP_Gf ) then
                 call GF_COP(r_oDev,Gf_tri,zwork_tri,pvt, &
-                     TSHS%sp,H,TSHS%sc_off, kpt, orb_J)
+                     TSHS%sp,H,TSHS%sc_off, kpt, phase, dev_M)
                 if ( dH%lvl > 0 ) then
-                   call GF_COHP_add_dH(dH%d, TSHS%sc_off, &
-                        kpt, Gf_tri, zwork_tri, r_oDev, orb_J, pvt)
+                   call GF_COHP_add_dH(dH%d, TSHS%sc_off, kpt, &
+                       phase, Gf_tri, zwork_tri, r_oDev, dev_M, pvt)
                 end if
-                call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COHP', orb_J)
+                call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COHP', dev_M)
              end if
 #endif
 
@@ -1033,7 +1064,7 @@
           ! ****************
           ! * Column Gf    *
           ! ****************
-          if ( 'T-Gf' .in. save_DATA ) then
+          if ( calc_T_Gf ) then
 
              ! We are allowed to calculate the transmission
              ! only by using the diagonal
@@ -1070,7 +1101,7 @@
                  call invert_BiasTriMat_rgn(GF_tri,zwork_tri, &
                       r_oDev, pvt, Elecs(iEl)%o_inD)
  
-                 if ( 'T-sum-out' .in. save_DATA ) then
+                 if ( calc_T_out ) then
                     call Gf_Gamma(zwork_tri,Elecs(iEl),T(N_Elec+1,iEl))
                  end if
  
@@ -1092,7 +1123,7 @@
               ! * calc A-matrix *
               ! *****************
               if ( .not. cE%fake ) then
-                 if ( 'T-sum-out' .in. save_DATA ) then
+                 if ( calc_T_out ) then
                     call dir_GF_Gamma_GF(Gf_tri, zwork_tri, r_oDev, pvt, &
                          Elecs(iEl), A_parts, &
                          TrGfG = T(N_Elec+1,iEl))
@@ -1114,20 +1145,25 @@
                 end if
                 
 #ifdef NCDF_4
+                if ( calc_DM_A ) then
+                   call A_DM(TSHS%sc_off,kpt,phase,zwork_tri,r_oDev,pvt, dev_M)
+                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'DM', dev_M, &
+                         Elecs(iEl))
+                end if
                 if ( calc_COOP_A ) then
                    call A_COP(r_oDev,zwork_tri,pvt, &
-                        TSHS%sp,S,TSHS%sc_off, kpt, orb_J)
-                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COOP', orb_J, &
+                        TSHS%sp,S,TSHS%sc_off, kpt, phase, dev_M)
+                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COOP', dev_M, &
                         Elecs(iEl))
                 end if
                 if ( calc_COHP_A ) then
                    call A_COP(r_oDev,zwork_tri,pvt, &
-                        TSHS%sp,H,TSHS%sc_off, kpt, orb_J)
+                        TSHS%sp,H,TSHS%sc_off, kpt, phase, dev_M)
                    if ( dH%lvl > 0 ) then
                       call A_COHP_add_dH(dH%d, TSHS%sc_off, &
-                           kpt, zwork_tri, r_oDev, orb_J, pvt)
+                           kpt, phase, zwork_tri, r_oDev, dev_M, pvt)
                    end if
-                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COHP', orb_J, &
+                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'COHP', dev_M, &
                         Elecs(iEl))
                 end if
                 
@@ -1135,21 +1171,21 @@
 
 #ifdef TBT_PHONON
                    call orb_current(TSHS%sp,H,S,TSHS%sc_off, &
-                        kpt, &
-                        cOmega,zwork_tri,r_oDev,orb_J,pvt)
+                        kpt, phase, &
+                        cOmega,zwork_tri,r_oDev,dev_M,pvt)
 #else
                    call orb_current(TSHS%sp,H,S,TSHS%sc_off, &
-                        kpt, &
-                        cE,zwork_tri,r_oDev,orb_J,pvt)
+                        kpt, phase, &
+                        cE,zwork_tri,r_oDev,dev_M,pvt)
 #endif
                    if ( dH%lvl > 0 ) then
                       call orb_current_add_dH(dH%d, TSHS%sc_off, &
-                           kpt, zwork_tri, r_oDev, orb_J, pvt)
+                           kpt, phase, zwork_tri, r_oDev, dev_M, pvt)
                    end if
 
                    ! We need to save it immediately, we
                    ! do not want to have several arrays in memory
-                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'J', orb_J, &
+                   call state_cdf_save_sp_dev(TBTcdf, ikpt, nE, 'J', dev_M, &
                         Elecs(iEl))
 
                 end if
@@ -1163,8 +1199,7 @@
                 ! calculate this.
                 if ( (.not. T_all) .and. &
                      jEl < iEl ) cycle
-                if ( ('T-sum-out' .nin. save_DATA ) .and. &
-                     iEl == jEl ) cycle
+                if ( (.not. calc_T_out) .and. iEl == jEl ) cycle
 
                 ! Notice that the Gf.G1.Gf.G2 can be performed
                 ! for all other electrodes as long as we
@@ -1210,7 +1245,7 @@
                save_DATA )
 #endif
 
-          call timer('DOS-Gf-A-T',2)
+          call timer('analysis',2)
 #ifdef NCDF_4
           end if ! .not. proj-only
           end if ! .not. Sigma-only
@@ -1221,7 +1256,7 @@
 
           if ( N_proj_T > 0 ) then
 
-          call timer('T-proj',1)
+          call timer('analysis-proj',1)
 
           ! Calculate the projections
           do ipt = 1 , N_proj_T
@@ -1236,7 +1271,7 @@
                 ! superfluous checks in the following block
                 if ( ('proj-orb-current' .in. save_DATA) .and. p_E%idx > 0 ) then
 
-                   call proj_cdf_save_J(PROJcdf, ikpt, nE, p_E, orb_J)
+                   call proj_cdf_save_J(PROJcdf, ikpt, nE, p_E, dev_M)
 
                 end if
 #endif
@@ -1305,7 +1340,7 @@
                         bTk(1+size(proj_T(ipt)%R),ipt))
                 end if
                else
-                if ( 'T-sum-out' .in. save_DATA ) then
+                if ( 'proj-T-sum-out' .in. save_DATA ) then
                    call dir_GF_Gamma_GF(Gf_tri, zwork_tri, r_oDev, pvt, &
                         Elecs(iEl), proj_parts, &
                         TrGfG = bTk(1+size(proj_T(ipt)%R),ipt))
@@ -1336,22 +1371,22 @@
 
 #ifdef TBT_PHONON
                   call orb_current(TSHS%sp,H,S,TSHS%sc_off, &
-                       kpt, &
-                       cOmega,zwork_tri,r_oDev,orb_J,pvt)
+                       kpt, phase, &
+                       cOmega,zwork_tri,r_oDev,dev_M,pvt)
 #else
                   call orb_current(TSHS%sp,H,S,TSHS%sc_off, &
-                       kpt, &
-                       cE,zwork_tri,r_oDev,orb_J,pvt)
+                       kpt, phase, &
+                       cE,zwork_tri,r_oDev,dev_M,pvt)
 #endif
                   if ( dH%lvl > 0 ) then
                      call orb_current_add_dH(dH%d, TSHS%sc_off, &
-                          kpt, zwork_tri, r_oDev, orb_J, pvt)
+                          kpt, phase, zwork_tri, r_oDev, dev_M, pvt)
                   end if
                      
                   ! We need to save it immediately, we
                   ! do not want to have several arrays in the
                   ! memory
-                  call proj_cdf_save_J(PROJcdf, ikpt, nE, p_E, orb_J)
+                  call proj_cdf_save_J(PROJcdf, ikpt, nE, p_E, dev_M)
 
                end if
 #endif
@@ -1426,7 +1461,7 @@
                ikpt,nE,N_proj_T,proj_T, &
                pDOS, bTk, N_eigen, bTkeig, save_DATA )
 
-          call timer('T-proj',2)
+          call timer('analysis-proj',2)
           
           end if
 #endif
@@ -1514,7 +1549,7 @@
        deallocate(bTkeig)
     end if
 
-    call delete(orb_J)
+    call delete(dev_M)
 
     ! Before we delete the Gf tri-diagonal matrix
     ! we need to create the sigma mean if requested.

=== modified file 'Util/TS/TBtrans/m_tbtrans.F90'
--- Util/TS/TBtrans/m_tbtrans.F90	2018-04-09 12:39:44 +0000
+++ Util/TS/TBtrans/m_tbtrans.F90	2018-05-04 19:08:45 +0000
@@ -301,18 +301,18 @@
           call name_save( ispin, TSHS%nspin,cdf_fname, end = 'nc')
           call init_cdf_save(cdf_fname,TSHS,r_oDev,DevTri,ispin, &
               N_Elec, Elecs, r_oElpd, ElTri, &
-              nkpt, kpt, wkpt, NEn, r_aDev, r_aBuf, sp_dev_sc, save_DATA )
+              nkpt, kpt, wkpt, NEn, tbt_Eta, r_aDev, r_aBuf, sp_dev_sc, save_DATA )
        end if
        
        call name_save( ispin, TSHS%nspin,cdf_fname_sigma, end = 'SE.nc')
        call init_Sigma_save(cdf_fname_sigma,TSHS,r_oDev,DevTri,ispin,N_Elec, Elecs, &
-            nkpt, kpt, wkpt, NEn, r_aDev, r_aBuf )
+            nkpt, kpt, wkpt, NEn, tbt_Eta, r_aDev, r_aBuf )
 
        if ( ('Sigma-only'.nin.save_DATA) ) then
        
           call name_save( ispin, TSHS%nspin, cdf_fname_proj, end = 'Proj.nc' )
           call init_Proj_save( cdf_fname_proj,TSHS,r_oDev,DevTri,ispin, N_Elec, Elecs, &
-               nkpt, kpt, wkpt, NEn , r_aDev, r_aBuf, sp_dev_sc, save_DATA )
+               nkpt, kpt, wkpt, NEn, tbt_Eta, r_aDev, r_aBuf, sp_dev_sc, save_DATA )
        end if
 
        if ( n_k /= 0 ) then

=== modified file 'version.info'
--- version.info	2018-04-27 09:35:52 +0000
+++ version.info	2018-05-04 19:08:45 +0000
@@ -1,2 +1,1 @@
-siesta-4.1--903
-
+siesta-4.1--903--ts-10