diff --git a/.gitignore b/.gitignore
index d1a19546..d42d53c7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,3 +11,4 @@ __pycache__
 _gitmsg.saved.txt
 *.egg-info
 dist
+.vscode/
diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt
index 03df19ba..f0e21a89 100644
--- a/slateratom/lib/CMakeLists.txt
+++ b/slateratom/lib/CMakeLists.txt
@@ -31,6 +31,7 @@ set(sources-fpp
   blasroutines.F90
   broydenmixer.F90
   confinement.F90
+  diismixer.F90
   input.F90
   lapackroutines.F90
   simplemixer.F90
diff --git a/slateratom/lib/broydenmixer.F90 b/slateratom/lib/broydenmixer.F90
index b29a7e03..4251df08 100644
--- a/slateratom/lib/broydenmixer.F90
+++ b/slateratom/lib/broydenmixer.F90
@@ -48,11 +48,11 @@ module broydenmixer
     !> Weights for prev. iterations
     real(dp), allocatable :: ww(:)
 
-    !> Charge difference in last iteration
-    real(dp), allocatable :: qDiffLast(:)
+    !> Input increment in last iteration
+    real(dp), allocatable :: incrementLast(:)
 
-    !> Input charge in last iteration
-    real(dp), allocatable :: qInpLast(:)
+    !> Input quantity that is being mixed, last iteration
+    real(dp), allocatable :: inputLast(:)
 
     !> Storage for the "a" matrix
     real(dp), allocatable :: aa(:,:)
@@ -70,7 +70,7 @@ module broydenmixer
 
   !> Creates a Broyden mixer instance.
   !! The weight associated with an iteration is calculated as weigthFac/ww where ww is the Euclidean
-  !! norm of the charge difference vector. If the calculated weigth is outside of the
+  !! norm of the quantity difference vector. If the calculated weigth is outside of the
   !! [minWeight, maxWeight] region it is replaced with the appropriate boundary value.
   subroutine TBroydenMixer_init(this, mIter, mixParam, omega0, minWeight, maxWeight, weightFac)
 
@@ -107,8 +107,8 @@ subroutine TBroydenMixer_init(this, mIter, mixParam, omega0, minWeight, maxWeigh
     this%maxWeight = maxWeight
     this%weightFac = weightFac
     allocate(this%ww(mIter-1))
-    allocate(this%qInpLast(this%nElem))
-    allocate(this%qDiffLast(this%nElem))
+    allocate(this%inputLast(this%nElem))
+    allocate(this%incrementLast(this%nElem))
     allocate(this%aa(mIter-1, mIter-1))
     allocate(this%dF(this%nElem, mIter - 1))
     allocate(this%uu(this%nElem, mIter - 1))
@@ -129,10 +129,10 @@ subroutine TBroydenMixer_reset(this, nElem)
 
     if (nElem /= this%nElem) then
       this%nElem = nElem
-      deallocate(this%qInpLast)
-      deallocate(this%qDiffLast)
-      allocate(this%qInpLast(this%nElem))
-      allocate(this%qDiffLast(this%nElem))
+      deallocate(this%inputLast)
+      deallocate(this%incrementLast)
+      allocate(this%inputLast(this%nElem))
+      allocate(this%incrementLast(this%nElem))
       deallocate(this%dF)
       allocate(this%dF(this%nElem, this%mIter - 1))
       deallocate(this%uu)
@@ -153,24 +153,24 @@ end subroutine TBroydenMixer_reset
   !!   The restriction arises from the assumption that the dot-products of density matrices are
   !!   real-valued (imaginary parts add up to zero due to the hermitian property) and the linear
   !!   system of equations remains real-valued.
-  subroutine TBroydenMixer_mix(this, qInpResult, qDiff)
+  subroutine TBroydenMixer_mix(this, inputResult, increment)
 
     !> The Broyden mixer
     type(TBroydenMixer), intent(inout) :: this
 
     !> Input charges on entry, mixed charges on exit
-    real(dp), intent(inout) :: qInpResult(:)
+    real(dp), intent(inout) :: inputResult(:)
 
     !> Charge difference between output and input charges
-    real(dp), intent(in) :: qDiff(:)
+    real(dp), intent(in) :: increment(:)
 
     this%iIter = this%iIter + 1
     if (this%iIter > this%mIter) then
       error stop "Broyden mixer: Maximal nr. of steps exceeded"
     end if
 
-    call modifiedBroydenMixing(qInpResult, this%qInpLast, this%qDiffLast, this%aa,&
-        & this%ww, this%iIter, qDiff, this%alpha, this%omega0, this%minWeight, this%maxWeight,&
+    call modifiedBroydenMixing(inputResult, this%inputLast, this%incrementLast, this%aa,&
+        & this%ww, this%iIter, increment, this%alpha, this%omega0, this%minWeight, this%maxWeight,&
         & this%weightFac, this%nElem, this%dF, this%uu)
 
   end subroutine TBroydenMixer_mix
diff --git a/slateratom/lib/diagonalizations.f90 b/slateratom/lib/diagonalizations.f90
index 2fd4a10a..9aae3d5b 100644
--- a/slateratom/lib/diagonalizations.f90
+++ b/slateratom/lib/diagonalizations.f90
@@ -14,7 +14,7 @@ module diagonalizations
 
   !> Diagonalizes overlap matrix to check for linear dependency of basis set.
   !! Implicitely LAPACK's dsyev is called.
-  subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
+  subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss, invsqrt_ss)
 
     !> maximum angular momentum
     integer, intent(in) :: max_l
@@ -28,6 +28,9 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
     !> overlap supervector
     real(dp), intent(in) :: ss(0:, :,:)
 
+    !> inv. sqrt. of the overlap supervector
+    real(dp), intent(inout) :: invsqrt_ss(0:, :,:)
+
     !! overlap matrices
     real(dp), allocatable :: overlap(:,:)
 
@@ -35,7 +38,7 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
     real(dp), allocatable :: eigenvalues(:)
 
     !> auxiliary variables
-    integer :: ll, diagsize
+    integer :: ll, diagsize, ii
 
     do ll = 0, max_l
 
@@ -46,7 +49,7 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
 
       overlap = ss(ll, :,:)
 
-      call heev(overlap, eigenvalues, 'U', 'N')
+      call heev(overlap, eigenvalues, 'U', 'V')
 
       write(*, '(A,I3,A,E16.8)') 'Smallest eigenvalue of overlap for l= ', ll, ' : ', eigenvalues(1)
 
@@ -57,6 +60,12 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
         stop
       end if
 
+      ! compute S^(-1/2)
+      do ii = 1, diagsize
+        invsqrt_ss(ll, ii, ii) = eigenvalues(ii)**(-1.0_dp / 2)
+      end do
+      invsqrt_ss(ll,:,:) = matmul(overlap, matmul(invsqrt_ss(ll,:,:), transpose(overlap)))
+
       deallocate(overlap, eigenvalues)
 
     end do
diff --git a/slateratom/lib/diismixer.F90 b/slateratom/lib/diismixer.F90
new file mode 100644
index 00000000..0765339c
--- /dev/null
+++ b/slateratom/lib/diismixer.F90
@@ -0,0 +1,219 @@
+#:include 'common.fypp'
+
+!> Contains the DIIS mixer
+!! The DIIS mixing is done by building a weighted combination over the previous input charges to
+!! minimise the residue of the error.
+!! Only a specified number of previous charge vectors are considered.
+!! The modification based on from Kovalenko et al. (J. Comput. Chem., 20: 928-936 1999) and Patrick
+!! Briddon to add a contribution from the gradient vector as well is also used.
+!! In order to use the mixer you have to create and reset it.
+!!
+!! Code is adapted from DFTB+.
+module diismixer
+  use common_accuracy, only : dp
+  use lapackroutines, only : gesv
+  implicit none
+
+  private
+
+  public :: TDiisMixer, TDiisMixer_init, TDiisMixer_mix, TDiisMixer_reset
+
+  !> Contains the necessary data for an DIIS mixer.
+  type TDiisMixer
+    private
+
+    !> Initial mixing parameter
+    real(dp) :: initMixParam
+
+    !> Max. nr. of stored prev. vectors
+    integer :: mPrevVector
+
+    !> Nr. of stored previous vectors
+    integer :: iPrevVector
+
+    !> Nr. of elements in the vectors
+    integer :: nElem
+
+    !> Index for the storage
+    integer :: indx
+
+    !> Stored previous input quantities
+    real(dp), allocatable :: prevInput(:,:)
+
+    !> Stored differences of previous input quantities
+    real(dp), allocatable :: prevIncrement(:,:)
+
+    !> Stored prev. error vectors
+    real(dp), allocatable :: prevErrorVec(:,:)
+
+    !> True if DIIS used from iteration 2 as well as mixing
+    logical :: tFromStart
+
+    !> Alpha factor to add in new information
+    real(dp) :: alpha
+
+    contains
+      procedure :: reset => TDiisMixer_reset
+      procedure :: mix1D => TDiisMixer_mix
+  end type TDiisMixer
+
+contains
+
+  !> Initializes a DIIS mixer instance.
+  subroutine TDiisMixer_init(this, iGenerations, initMixParam, tFromStart)
+
+    !> Pointer to an initialized DIIS mixer on exit
+    type(TDiisMixer), intent(out) :: this
+
+    !> Number of generations to consider (including current)
+    integer, intent(in) :: iGenerations
+
+    !> Damping parameter for the first mixing steps
+    real(dp), intent(in) :: initMixParam
+
+    !> Use DIIS from step 2 onwards?
+    logical, intent(in) :: tFromStart
+
+    @:ASSERT(iGenerations >= 2)
+
+    this%nElem = 0
+    this%mPrevVector = iGenerations
+
+    allocate(this%prevInput(this%nElem, this%mPrevVector))
+    allocate(this%prevErrorVec(this%nElem, this%mPrevVector))
+    allocate(this%prevIncrement(this%nElem, this%mPrevVector))
+
+    this%initMixParam = initMixParam
+    this%tFromStart = tFromStart
+
+  end subroutine TDiisMixer_init
+
+
+  !> Makes the mixer ready for a new SCC cycle.
+  subroutine TDiisMixer_reset(this, nElem)
+
+    !> DIIS mixer instance
+    class(TDiisMixer), intent(inout) :: this
+
+    !> Nr. of elements in the vectors to mix
+    integer, intent(in) :: nElem
+
+    @:ASSERT(nElem > 0)
+
+    if (nElem /= this%nElem) then
+      this%nElem = nElem
+      deallocate(this%prevInput)
+      deallocate(this%prevErrorVec)
+      deallocate(this%prevIncrement)
+      allocate(this%prevInput(this%nElem, this%mPrevVector))
+      allocate(this%prevErrorVec(this%nElem, this%mPrevVector))
+      allocate(this%prevIncrement(this%nElem, this%mPrevVector))
+    end if
+    this%iPrevVector = 0
+    this%indx = 0
+
+  end subroutine TDiisMixer_reset
+
+
+  !> Mixes quantities according to the DIIS method.
+  subroutine TDiisMixer_mix(this, inputResult, inputIncrement, errorVector)
+
+    !> Pointer to the diis mixer
+    class(TDiisMixer), intent(inout) :: this
+
+    !> Input quantity on entry, mixed quantity on exit.
+    real(dp), intent(inout) :: inputResult(:)
+
+    !> Increment vector between input and output quantities
+    real(dp), intent(in) :: inputIncrement(:)
+
+    !> Error metric vector
+    real(dp), intent(in) :: errorVector(:)
+
+
+    real(dp), allocatable :: aa(:,:), bb(:,:)
+    integer :: ii, jj
+
+    @:ASSERT(size(inputResult) == this%nElem)
+    @:ASSERT(size(errorVector) == this%nElem)
+
+    if (this%iPrevVector < this%mPrevVector) then
+      this%iPrevVector = this%iPrevVector + 1
+    end if
+
+    call storeVectors(this%prevInput, this%prevIncrement, this%prevErrorVec, this%indx,&
+        & inputResult, inputIncrement, errorVector, this%mPrevVector)
+    if (this%tFromStart .or. this%iPrevVector == this%mPrevVector) then
+
+      allocate(aa(this%iPrevVector + 1, this%iPrevVector + 1))
+      allocate(bb(this%iPrevVector + 1, 1))
+
+      aa(:,:) = 0.0_dp
+      bb(:,:) = 0.0_dp
+
+      ! (due to the hermitian property of our density matrices, the dot-product below is real)
+      do ii = 1, this%iPrevVector
+        do jj = 1, this%iPrevVector
+          aa(ii, jj) = dot_product(this%prevErrorVec(:, ii), this%prevErrorVec(:, jj))
+        end do
+      end do
+      aa(this%iPrevVector + 1, 1:this%iPrevVector) = -1.0_dp
+      aa(1:this%iPrevVector, this%iPrevVector + 1) = -1.0_dp
+
+      bb(this%iPrevVector + 1, 1) = -1.0_dp
+
+      ! Solve DIIS system of linear equations
+      call gesv(aa, bb)
+
+      inputResult(:) = 0.0_dp
+      do ii = 1, this%iPrevVector
+        inputResult(:) = inputResult + bb(ii, 1) * (this%prevInput(:, ii) + this%prevIncrement(:, ii))
+      end do
+
+    end if
+
+    if (this%iPrevVector < this%mPrevVector) then
+      ! First few iterations return simple mixed vector
+      inputResult(:) = inputResult + this%initMixParam * inputIncrement(:)
+    end if
+
+  end subroutine TDiisMixer_mix
+
+
+  !> Stores a vector pair in a limited storage.
+  !! If the stack is full, oldest vector pair is overwritten.
+  subroutine storeVectors(prevInp, prevIncrement, prevErrorVector, indx, input, increment,&
+        & errorVector, mPrevVector)
+
+    !> Contains previous vectors of the first type
+    real(dp), intent(inout) :: prevInp(:,:)
+
+    !> Contains previous vectors of the second type
+    real(dp), intent(inout) :: prevErrorVector(:,:)
+
+    !> Contains previous differences of vectors of first type
+    real(dp), intent(inout) :: prevIncrement(:,:)
+
+    !> Indexing of data
+    integer, intent(inout) :: indx
+
+    !> New first vector
+    real(dp), intent(in) :: input(:)
+
+    !> New increment of first vector
+    real(dp), intent(in) :: increment(:)
+
+    !> New second vector
+    real(dp), intent(in) :: errorVector(:)
+
+    !> Size of the stacks.
+    integer, intent(in) :: mPrevVector
+
+    indx = mod(indx, mPrevVector) + 1
+    prevInp(:, indx) = input
+    prevIncrement(:, indx) = increment
+    prevErrorVector(:, indx) = errorVector
+    
+  end subroutine storeVectors
+
+end module diismixer
\ No newline at end of file
diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90
index 92e50a35..8da01e17 100644
--- a/slateratom/lib/globals.f90
+++ b/slateratom/lib/globals.f90
@@ -5,6 +5,7 @@ module globals
   use mixer, only : TMixer, TMixer_init, TMixer_reset, mixerTypes
   use broydenmixer, only : TBroydenMixer, TBroydenMixer_init
   use simplemixer, only : TSimpleMixer, TSimpleMixer_init
+  use diismixer, only : TDiisMixer, TDiisMixer_init
   use confinement, only : TConfInp
 
   implicit none
@@ -61,6 +62,9 @@ module globals
   !> overlap supervector
   real(dp), allocatable :: ss(:,:,:)
 
+  !> inv. square root of the overlap supervector
+  real(dp), allocatable :: invsqrt_ss(:,:,:)
+
   !> nucleus-electron supervector
   real(dp), allocatable :: uu(:,:,:)
 
@@ -85,8 +89,14 @@ module globals
   !> wavefunction coefficients
   real(dp), allocatable :: cof(:,:,:,:)
 
-  !> relative changes during scf
-  real(dp) :: change_max
+  !> max abs. change in potential matrix in the last scf step
+  real(dp) :: d_pot_max
+
+  !> max abs. change in occupied Fock eigenvalues in the last scf step
+  real(dp) :: d_spectrum_max
+
+  !> max abs. value in the commutator from the current scf step
+  real(dp) :: commutator_max
 
   !> density matrix supervector
   real(dp), allocatable :: pp(:,:,:,:)
@@ -94,15 +104,21 @@ module globals
   !> fock matrix supervector
   real(dp), allocatable :: ff(:,:,:,:)
 
+  !> commutator [F, PS]
+  real(dp), allocatable :: commutator(:,:,:,:)
+
   !> potential matrix supervectors
   real(dp), allocatable :: pot_new(:,:,:,:), pot_old(:,:,:,:)
 
   !> eigenvalues
-  real(dp), allocatable :: eigval(:,:,:)
+  real(dp), allocatable :: eigval(:,:,:), eigval_old(:,:,:)
 
   !> zora scaled eigenvalues
   real(dp), allocatable :: eigval_scaled(:,:,:)
 
+  !> true, if SCF cycle reached convergency on a given quantity
+  logical :: tCommutatorConverged, tEnergyConverged, tSpectrumConverged
+
   !> total energy
   real(dp) :: total_ene
 
@@ -187,6 +203,9 @@ module globals
   !> broyden mixer (if used)
   type(TBroydenMixer), allocatable :: pBroydenMixer
 
+  !> DIIS mixer (if used)
+  type(TDiisMixer), allocatable :: pDiisMixer
+
   !> mixing factor
   real(dp) :: mixing_factor
 
@@ -220,16 +239,19 @@ subroutine allocate_globals()
 
     allocate(ss(0:max_l, problemsize, problemsize))
     write(*, '(A,I0,A)') 'Size of one Supervectors is ', size(ss), ' double precision elements'
+    allocate(invsqrt_ss(0:max_l, problemsize, problemsize))
 
     allocate(uu(0:max_l, problemsize, problemsize))
     allocate(tt(0:max_l, problemsize, problemsize))
     allocate(vconf(num_mesh_points, 0:max_l))
     allocate(vconf_matrix(0:max_l, problemsize, problemsize))
     allocate(ff(2, 0:max_l, problemsize, problemsize))
+    allocate(commutator(2, 0:max_l, problemsize, problemsize))
     allocate(pot_old(2, 0:max_l, problemsize, problemsize))
     allocate(pot_new(2, 0:max_l, problemsize, problemsize))
 
     allocate(eigval(2, 0:max_l, problemsize))
+    allocate(eigval_old(2, 0:max_l, problemsize))
     allocate(eigval_scaled(2, 0:max_l, problemsize))
 
     allocate(jj(0:max_l, problemsize, problemsize, 0:max_l, problemsize, problemsize))
@@ -253,6 +275,7 @@ subroutine allocate_globals()
     ddrho(:,:) = 0.0_dp
 
     eigval(:,:,:) = 0.0_dp
+    eigval_old(:,:,:) = 0.0_dp
     eigval_scaled(:,:,:) = 0.0_dp
 
     cof(:,:,:,:) = 0.0_dp
@@ -273,6 +296,10 @@ subroutine allocate_globals()
       call TBroydenMixer_init(pBroydenMixer, maxiter, mixing_factor, 0.01_dp, 1.0_dp, 1.0e5_dp,&
           & 1.0e-2_dp)
       call TMixer_init(pMixer, pBroydenMixer)
+    case(mixerTypes%diis)
+      allocate(pDiisMixer)
+      call TDiisMixer_init(pDiisMixer, 10, mixing_factor, .false.)
+      call TMixer_init(pMixer, pDiisMixer)
     case default
       error stop "Unknown mixer type."
     end select
diff --git a/slateratom/lib/hamiltonian.f90 b/slateratom/lib/hamiltonian.f90
index dab8d521..5d354419 100644
--- a/slateratom/lib/hamiltonian.f90
+++ b/slateratom/lib/hamiltonian.f90
@@ -3,9 +3,10 @@ module hamiltonian
 
   use common_accuracy, only : dp
   use dft, only : dft_exc_matrixelement
-  use mixer, only : TMixer, TMixer_mix
-  use zora_routines, only : zora_t_correction
+  use mixer, only : TMixer, TMixer_mix, TMixer_getMixerType, TMixer_reset
+  use utilities, only : compute_commutator
   use xcfunctionals, only : xcFunctional
+  use zora_routines, only : zora_t_correction
 
   implicit none
   private
@@ -19,7 +20,7 @@ module hamiltonian
   !> Main driver routine for Fock matrix build-up. Also calls mixer with potential matrix.
   subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l,&
       & num_alpha, poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha,&
-      & pot_old, pot_new, tZora, ff, camAlpha, camBeta)
+      & pot_old, pot_new, overlap, invsqrt_ss, commutator, tZora, ff, camAlpha, camBeta)
 
     !> mixer instances
     type(TMixer), intent(inout) :: pMixer
@@ -87,6 +88,15 @@ subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp
     !> new potential
     real(dp), intent(out) :: pot_new(:,0:,:,:)
 
+    !> Overlap matrix S
+    real(dp), intent(in) :: overlap(0:, :,:)
+    
+    !> S^(-1/2)
+    real(dp), intent(in) :: invsqrt_ss(0:, :,:)
+
+    !> commutator S^(-1/2) [F,PS] S^(-1/2)
+    real(dp), intent(out) :: commutator(:, 0:, :, :)
+
     !> true, if zero-order regular approximation for relativistic effects is desired
     logical, intent(in) :: tZora
 
@@ -179,10 +189,52 @@ subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp
     pot_new(1, :,:,:) = - real(nuc, dp) * uu + j_matrix - k_matrix(1, :,:,:)
     pot_new(2, :,:,:) = - real(nuc, dp) * uu + j_matrix - k_matrix(2, :,:,:)
 
+    ! pre-build Fock matrix for commutator computation
+    do ii = 0, max_l
+      ss = 0
+      do jjj = 1, num_alpha(ii)
+        do kkk = 1, poly_order(ii)
+          ss = ss + 1
+          ttt = 0
+          do ll = 1, num_alpha(ii)
+            do mm = 1, poly_order(ii)
+              ttt = ttt + 1
+
+              ff(1, ii, ss, ttt) = tt(ii, ss, ttt) + pot_new(1, ii, ss, ttt) + vconf(ii, ss, ttt)
+              ff(2, ii, ss, ttt) = tt(ii, ss, ttt) + pot_new(2, ii, ss, ttt) + vconf(ii, ss, ttt)
+
+              if (tZora) then
+                ff(1, ii, ss, ttt) = ff(1, ii, ss, ttt) + t_zora(1, ii, ss, ttt)
+                ff(2, ii, ss, ttt) = ff(2, ii, ss, ttt) + t_zora(2, ii, ss, ttt)
+              end if
+
+            end do
+          end do
+        end do
+      end do
+    end do
+
+    ! compute S^(-1/2) [F,PS] S^(-1/2)
+    call compute_commutator(max_l, num_alpha, poly_order, ff, pp, overlap, invsqrt_ss, commutator)
+
+    ! Not sure: before or after mixer (potential .ne. Matrix elements)?
+    ! Should be irrelevant once self-consistency is reached.
+    if (tZora .and. (iScf /= 0)) then
+      call zora_t_correction(1, t_zora, max_l, num_alpha, alpha, poly_order, num_mesh_points,&
+          & weight, abcissa, vxc, nuc, pp, problemsize)
+    end if
+
     ! mixer
     allocate(pot_diff, mold=pot_old)
-    pot_diff(:,0:,:,:) = pot_old - pot_new
-    call TMixer_mix(pMixer, pot_new, pot_diff)
+    ! pot_diff(:,0:,:,:) = pot_old - pot_new
+    pot_diff(:,0:,:,:) = pot_new - pot_old
+
+    call TMixer_mix(pMixer, pot_new, pot_diff, commutator)
+
+    ! guard against uninitalised arrays on step 0
+    if (iScf == 0) then
+      call TMixer_reset(pMixer, size(pot_new))
+    end if
 
     ! Not sure: before or after mixer (potential .ne. Matrix elements)?
     ! Should be irrelevant once self-consistency is reached.
@@ -216,6 +268,8 @@ subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp
       end do
     end do
 
+   deallocate(pot_diff)
+
   end subroutine build_hamiltonian
 
 
diff --git a/slateratom/lib/input.F90 b/slateratom/lib/input.F90
index c3f01454..d26677a3 100644
--- a/slateratom/lib/input.F90
+++ b/slateratom/lib/input.F90
@@ -252,7 +252,7 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     read(*,*) tPrintEigvecs
 
     write(*, '(A)') 'Enter mixer and mixing parameter <1:&
-        & 1: Simple mixer, 2: Broyden mixer'
+        & 1: Simple mixer, 2: Broyden mixer, 3: DIIS'
     read(*,*) mixnr, mixing_factor
 
   end subroutine read_input_1
diff --git a/slateratom/lib/lapackroutines.F90 b/slateratom/lib/lapackroutines.F90
index 390ba7cc..fe8836d0 100644
--- a/slateratom/lib/lapackroutines.F90
+++ b/slateratom/lib/lapackroutines.F90
@@ -4,11 +4,12 @@
 !! The interface of all LAPACK calls must be defined in the module lapack.
 module lapackroutines
 
-  use common_accuracy, only : dp, rdp
+  use common_accuracy, only : rdp
+  use common_message, only : error
   implicit none
 
   private
-  public :: getrf, getrs
+  public :: gesv, getrf, getrs
 
 
   !> Computes the LU decomposition of a general rectangular matrix using partial pivoting with row
@@ -29,8 +30,70 @@ module lapackroutines
   end interface getrs
 
 
+  !> Computes the solution to a real system of linear equations A * X = B, where A is an N-by-N
+  !> matrix and X and B are N-by-NRHS matrices
+  interface gesv
+    module procedure gesv_dble
+  end interface gesv
+
 contains
 
+  !> Solves a general system of linear equations A * X = B, where A is a real matrix.
+  subroutine gesv_dble(aa, bb, nEquation, nSolution)
+
+    !> Contains the coefficients on entry, the LU factorisation on exit.
+    real(rdp), intent(inout) :: aa(:,:)
+
+    !> Right hand side(s) of the linear equation on entry, solution(s) on exit.
+    real(rdp), intent(inout) :: bb(:,:)
+
+    !> The size of the problem (nr. of variables and equations). Must be only specified if different
+    !> from size(aa, dim=1).
+    integer, intent(in), optional :: nEquation
+
+    !> Nr. of right hand sides (nr. of solutions). Must be only specified if different from size(b,
+    !> dim=2).
+    integer, intent(in), optional :: nSolution
+
+    integer :: info
+    integer :: nn, nrhs, lda, ldb
+    integer, allocatable :: ipiv(:)
+    character(len=100) :: error_string
+
+    lda = size(aa, dim=1)
+    if (present(nEquation)) then
+      @:ASSERT(nEquation >= 1 .and. nEquation <= lda)
+      nn = nEquation
+    else
+      nn = lda
+    end if
+    @:ASSERT(size(aa, dim=2) >= nn)
+
+    ldb = size(bb, dim=1)
+    @:ASSERT(ldb >= nn)
+    nrhs = size(bb, dim=2)
+    if (present(nSolution)) then
+      @:ASSERT(nSolution <= nrhs)
+      nrhs = nSolution
+    end if
+
+    info = 0
+    allocate(ipiv(nn))
+    call dgesv(nn, nrhs, aa, lda, ipiv, bb, ldb, info)
+
+    if (info /= 0) then
+      if (info < 0) then
+        write(error_string, "(A,I0)")'Failure in dgesv illegal argument at position : ', info
+      else
+        write(error_string, "(A,I0)")'Linear dependent system in dgesv,&
+            & info flag : ', info
+      end if
+      call error(error_string)
+    end if
+
+  end subroutine gesv_dble
+
+
   !> Double precision version of getrf.
   subroutine getrf_dble(aa, ipiv, nRow, nColumn, iError)
 
@@ -162,4 +225,30 @@ subroutine getrs1_dble(amat, ipiv, bvec, trans, iError)
 
   end subroutine getrs1_dble
 
+  !> Solves overdetermined or underdetermined systems for GE matrices
+  ! subroutine dgelss_dble(amat, bmat, trans)
+
+  !   !> Matrix of the linear system
+  !   real(rdp), intent(inout) :: amat(:,:)
+
+  !   !>
+  !   real(rdp), intent(inout) :: bmat
+
+  !   !> Optional transpose (defaults to 'n')
+  !   character(len=1), intent(in), optional :: trans
+
+  !   !> Error flag, zero on successful exit
+  !   integer, intent(out), optional :: iError
+
+  !   integer :: nn, lda, info, lwork
+  !   real(rdp), allocatable :: work(:)
+  !   real(rdp) :: work2(1)
+    
+  !   integer :: info
+  !   allocate(svals(n), work(8*n))
+  !   call dgelss(n, n, 1, amat, size(amat,1), bmat, size(bmat,1), svals, 1.0e-10_dp, rank, work, 8*n, info)
+  !   if (info /= 0) call error("DIIS: dgelss failed to converge")
+  !     deallocate(svals, work)
+  ! end subroutine dgelss_dble
+
 end module lapackroutines
diff --git a/slateratom/lib/mixer.f90 b/slateratom/lib/mixer.f90
index 1ac8230f..35f7ab61 100644
--- a/slateratom/lib/mixer.f90
+++ b/slateratom/lib/mixer.f90
@@ -4,17 +4,18 @@ module mixer
   use common_accuracy, only : dp
   use broydenmixer, only : TBroydenMixer, TBroydenMixer_mix, TBroydenMixer_reset
   use simplemixer, only : TSimpleMixer, TSimpleMixer_mix, TSimpleMixer_reset
+  use diismixer, only : TDiisMixer, TDiisMixer_init, TDiisMixer_mix, TDiisMixer_reset
   implicit none
 
   private
-  public :: TMixer, TMixer_init, TMixer_reset, TMixer_mix, mixerTypes
+  public :: TMixer, TMixer_init, TMixer_getMixerType, TMixer_reset, TMixer_mix, mixerTypes
 
 
   !> Interface type for mixers
   type TMixer
     private
 
-    !> Numerical type of mixer 1:2
+    !> Numerical type of mixer 1:3
     integer :: mixerType
 
     !> Simple mixer instance
@@ -23,6 +24,9 @@ module mixer
     !> Broyden mixer instance
     type(TBroydenMixer), allocatable :: pBroydenMixer
 
+    !> DIIS mixer instance
+    type(TDiisMixer), allocatable :: pDiisMixer
+
   end type TMixer
 
 
@@ -30,6 +34,7 @@ module mixer
   interface TMixer_init
     module procedure TMixer_initSimple
     module procedure TMixer_initBroyden
+    module procedure TMixer_initDiis
   end interface TMixer_init
 
 
@@ -43,6 +48,7 @@ module mixer
   type :: TMixerTypesEnum
     integer :: simple = 1
     integer :: broyden = 2
+    integer :: diis = 3
   end type TMixerTypesEnum
 
   !> Contains mixer types
@@ -81,6 +87,33 @@ subroutine TMixer_initBroyden(this, pBroyden)
   end subroutine TMixer_initBroyden
 
 
+  !> Initializes a Broyden mixer.
+  subroutine TMixer_initDiis(this, pDiis)
+
+    !> Mixer instance
+    type(TMixer), intent(out) :: this
+
+    !> A valid Broyden mixer instance on exit
+    type(TDiisMixer), allocatable, intent(inout) :: pDiis
+
+    this%mixerType = mixerTypes%diis
+    call move_alloc(pDiis, this%pDiisMixer)
+
+  end subroutine TMixer_initDiis
+
+
+  !> Returns mixer type as an int (1: simple, 2: Broyden)
+  pure function TMixer_getMixerType(this) result(res)
+    !> Mixer instance
+    type(TMixer), intent(in) :: this
+
+    integer :: res
+
+    res = this%mixerType
+
+  end function
+
+
   !> Resets the mixer.
   subroutine TMixer_reset(this, nElem)
 
@@ -95,13 +128,15 @@ subroutine TMixer_reset(this, nElem)
       call TSimpleMixer_reset(this%pSimpleMixer, nElem)
     case(mixerTypes%broyden)
       call TBroydenMixer_reset(this%pBroydenMixer, nElem)
+    case(mixerTypes%diis)
+      call TDiisMixer_reset(this%pDiisMixer, nElem)
     end select
 
   end subroutine TMixer_reset
 
 
   !> Mixes two vectors.
-  subroutine TMixer_mix1D(this, inp, diff)
+  subroutine TMixer_mix1D(this, inp, diff, errvec)
 
     !> Mixer instance
     type(TMixer), intent(inout) :: this
@@ -112,18 +147,23 @@ subroutine TMixer_mix1D(this, inp, diff)
     !> Difference between input and output vectors (measure of lack of convergence)
     real(dp), intent(in) :: diff(:)
 
+    !> Error vector: measure of lack of convergence
+    real(dp), intent(in) :: errvec(:)
+
     select case (this%mixerType)
     case(mixerTypes%simple)
       call TSimpleMixer_mix(this%pSimpleMixer, inp, diff)
     case(mixerTypes%broyden)
       call TBroydenMixer_mix(this%pBroydenMixer, inp, diff)
+    case(mixerTypes%diis)
+      call TDiisMixer_mix(this%pDiisMixer, inp, diff, errvec)
     end select
 
   end subroutine TMixer_mix1D
 
 
   !> Mixes two 4D matrices.
-  subroutine TMixer_mix4D(this, inp, diff)
+  subroutine TMixer_mix4D(this, inp, diff, errvec)
 
     !> Mixer instance
     type(TMixer), intent(inout) :: this
@@ -131,19 +171,26 @@ subroutine TMixer_mix4D(this, inp, diff)
     !> Input vector on entry, result vector on exit
     real(dp), intent(inout), contiguous, target :: inp(:,0:,:,:)
 
-    !> Difference between input and output vectors (measure of lack of convergence)
+    !> Difference between input and output vectors
     real(dp), intent(in), contiguous, target :: diff(:,0:,:,:)
 
+    !> Error vector: measure of lack of convergence
+    real(dp), intent(in), contiguous, target :: errvec(:,0:,:,:)
+
     !! Difference between input and output vectors (1D pointer)
     real(dp), pointer :: pDiff(:)
 
     !! Input vector on entry, result vector on exit (1D pointer)
     real(dp), pointer :: pInp(:)
 
+    !> Error vector: measure of lack of convergence (1D pointers)
+    real(dp), pointer :: pErrVec(:)
+
     pInp(1:size(inp)) => inp
     pDiff(1:size(diff)) => diff
+    pErrVec(1:size(diff)) => errvec
 
-    call TMixer_mix1D(this, pInp, pDiff)
+    call TMixer_mix1D(this, pInp, pDiff, pErrVec)
 
   end subroutine TMixer_mix4D
 
diff --git a/slateratom/lib/utilities.f90 b/slateratom/lib/utilities.f90
index 1b1f8a66..c8c4d905 100644
--- a/slateratom/lib/utilities.f90
+++ b/slateratom/lib/utilities.f90
@@ -6,7 +6,8 @@ module utilities
   implicit none
   private
 
-  public :: check_convergence, check_electron_number
+  public :: check_convergence_pot, check_electron_number, check_convergence_commutator,&
+      & check_convergence_spectrum, compute_commutator
   public :: vector_length, fak, zeroOutCpotOfEmptyDensitySpinChannels
 
 
@@ -32,8 +33,147 @@ pure subroutine zeroOutCpotOfEmptyDensitySpinChannels(rho, vc)
   end subroutine zeroOutCpotOfEmptyDensitySpinChannels
 
 
+  !> Compute the commutator [F,PS] in MO basis  
+  subroutine compute_commutator(max_l, num_alpha, poly_order, ff, pp, ss, invsqrt_ss, commutator)
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of exponents in each shell
+    integer, intent(in) :: num_alpha(0:)
+
+    !> highest polynomial order + l in each shell
+    integer, intent(in) :: poly_order(0:)
+
+    !> Fock matrix
+    real(dp), intent(in) :: ff(:, 0:, :, :)
+
+    !> density matrix supervector
+    real(dp), intent(in) :: pp(:,0:,:,:)
+
+    !> overlap supervector
+    real(dp), intent(in) :: ss(0:,:,:)
+
+    !> inv. sqrt. of the overlap supervector
+    real(dp), intent(in) :: invsqrt_ss(0:,:,:)
+
+    !> commutator [F,PS]
+    real(dp), intent(out) :: commutator(:, 0:, :, :)
+
+    !! auxilliary variables
+    integer :: iSpin, ll
+
+    commutator(:, :, :, :) = 0.0_dp
+
+    ! Compute [F,PS] = FPS - SPF for each spin & ang. mom. block,
+    ! then orthogonalise like S^(-1/2) [F,PS] S^(-1/2)
+    do iSpin = 1, 2
+      do ll = 0, max_l
+        commutator(iSpin, ll, :, :) = &
+        & matmul(ff(iSpin, ll, :, :), matmul(pp(iSpin, ll, :, :), ss(ll, :, :)))&
+        & - matmul(matmul(ss(ll, :, :), pp(iSpin, ll, :, :)), ff(iSpin, ll, :, :))
+
+        commutator(iSpin, ll, :, :) = &
+        & matmul(invsqrt_ss(ll, :, :), matmul(commutator(iSpin, ll, :, :), invsqrt_ss(ll, :, :)))
+      end do
+    end do
+
+  end subroutine 
+
+
+  !> Check convergence by finding the maximum absolute value of the commutator [F,PS]
+  pure subroutine check_convergence_commutator(commutator, scftol, iScf, commutator_max,&
+      & tConverged)
+
+    !> commutator [F, PS]
+    real(dp), intent(in) :: commutator(:, 0:, :, :)
+
+    !> scf tolerance, i.e. convergence criteria
+    real(dp), intent(in) :: scftol
+
+    !> current SCF step
+    integer, intent(in) :: iScf
+
+    !> orbital gradient norm value
+    real(dp), intent(out) :: commutator_max
+
+    !> true, if SCF converged
+    logical, intent(out) :: tConverged
+
+    commutator_max = maxval(abs(commutator))
+
+    tConverged = commutator_max < scftol
+
+    if (iScf < 3) then
+      tConverged = .false.
+    end if
+
+  end subroutine check_convergence_commutator
+
+
+  !> Checks convergence by evaluating change in the occupied part of the eigenspectrum
+  pure subroutine check_convergence_spectrum(max_l, num_alpha, poly_order, eigval_new,&
+      & eigval_old, occ, scftol, iScf, res, tConverged)
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of exponents in each shell
+    integer, intent(in) :: num_alpha(0:)
+  
+    !> highest polynomial order + l in each shell
+    integer, intent(in) :: poly_order(0:)
+
+    !> old and new eigenspectra to compare
+    real(dp), intent(in) :: eigval_new(:,0:,:), eigval_old(:,0:,:)
+
+    !> occupations
+    real(dp), intent(in) :: occ(:,0:,:)
+
+    !> scf tolerance, i.e. convergence criteria
+    real(dp), intent(in) :: scftol
+
+    !> current SCF step
+    integer, intent(in) :: iScf
+
+    !> obtained change
+    real(dp), intent(out) :: res
+
+    !> true, if SCF converged
+    logical, intent(out) :: tConverged
+
+    ! Loop indices
+    integer :: iSpin, ll, diagsize, ii
+
+    ! Difference
+    real(dp) :: diff
+
+    ! Max. difference, only occupied orbitals contribute
+    res = 0.0_dp
+    do iSpin = 1, 2
+      do ll = 0, max_l
+        diagsize = num_alpha(ll) * poly_order(ll)
+        do ii = 1, diagsize
+          if (occ(iSpin, ll, ii) > 1e-16) then
+            diff = abs(eigval_new(iSpin, ll, ii) - eigval_old(iSpin, ll, ii)) 
+            if (diff > res) then
+              res = diff
+            end if
+          end if
+        end do
+      end do
+    end do
+
+    tConverged = res < scftol
+
+    if (iScf < 3) then
+      tConverged = .false.
+    end if
+
+  end subroutine check_convergence_spectrum
+
+
   !> Checks SCF convergence by comparing new and old potential.
-  pure subroutine check_convergence(pot_old, pot_new, max_l, problemsize, scftol, iScf, change_max,&
+  pure subroutine check_convergence_pot(pot_old, pot_new, max_l, problemsize, scftol, iScf, change_max,&
       & tConverged)
 
     !> old and new potential to compare
@@ -82,7 +222,7 @@ pure subroutine check_convergence(pot_old, pot_new, max_l, problemsize, scftol,
       tConverged = .false.
     end if
 
-  end subroutine check_convergence
+  end subroutine check_convergence_pot
 
 
   !> Checks conservation of electron number during SCF. If this fluctuates you are in deep trouble.
diff --git a/slateratom/prog/main.F90 b/slateratom/prog/main.F90
index 9ac24388..c8c5bc30 100644
--- a/slateratom/prog/main.F90
+++ b/slateratom/prog/main.F90
@@ -17,7 +17,8 @@ program HFAtom
       & write_waves_file_standard, write_wave_coeffs_file, cusp_values, writeAveragePotential
   use totalenergy, only : getTotalEnergy, getTotalEnergyZora
   use dft, only : check_accuracy, dft_start_pot, density_grid
-  use utilities, only : check_electron_number, check_convergence
+  use utilities, only : check_electron_number, check_convergence_pot, check_convergence_spectrum,&
+      & compute_commutator
   use zora_routines, only : scaled_zora
   use cmdargs, only : parse_command_arguments
   use common_poisson, only : TBeckeGridParams
@@ -94,6 +95,7 @@ program HFAtom
   ! build supervectors
   write(*, '(A)') 'Startup: Building Supervectors'
   call overlap(ss, max_l, num_alpha, alpha, poly_order)
+  
   call nuclear(uu, max_l, num_alpha, alpha, poly_order)
   call kinetic(tt, max_l, num_alpha, alpha, poly_order)
 
@@ -113,7 +115,7 @@ program HFAtom
   end if
 
   ! test for linear dependency
-  call diagonalize_overlap(max_l, num_alpha, poly_order, ss)
+  call diagonalize_overlap(max_l, num_alpha, poly_order, ss, invsqrt_ss)
 
   ! build supermatrices
   write(*, '(A)') 'Startup: Building Supermatrices'
@@ -147,15 +149,16 @@ program HFAtom
   ! kinetic energy, nuclear-electron, and confinement matrix elements which are constant during SCF
   call build_hamiltonian(pMixer, 0, tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
       & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,&
-      & pot_new, tZora, ff, camAlpha, camBeta)
+      & pot_new, ss, invsqrt_ss, commutator, tZora, ff, camAlpha, camBeta)
 
   ! self-consistency cycles
   write(*,*) 'Energies in Hartree'
   write(*,*)
-  write(*,*) ' Iter |   Total energy  |   HF-X energy  |   XC energy   |   Change in pot'
+  write(*,*) ' Iter |   Total energy  |  d(Spectrum)  |   d(Commutator)   |  d(Potential)'
   write(*,*) '--------------------------------------------------------------------------'
   lpScf: do iScf = 1, maxiter
 
+    eigval_old(:,:,:) = eigval
     pot_old(:,:,:,:) = pot_new
 
     ! diagonalize
@@ -171,7 +174,7 @@ program HFAtom
     ! build Fock matrix and get total energy during SCF
     call build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l,&
         & num_alpha, poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha,&
-        & pot_old, pot_new, tZora, ff, camAlpha, camBeta)
+        & pot_old, pot_new, ss, invsqrt_ss, commutator, tZora, ff, camAlpha, camBeta)
 
     if (tZora) then
       call getTotalEnergyZora(tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
@@ -185,10 +188,19 @@ program HFAtom
           & total_ene)
     end if
 
-    call check_convergence(pot_old, pot_new, max_l, problemsize, scftol, iScf, change_max,&
+    call check_convergence_pot(pot_old, pot_new, max_l, problemsize, scftol, iScf, d_pot_max,&
         & tConverged)
+    call check_convergence_spectrum(max_l, num_alpha, poly_order, eigval, eigval_old, occ,&
+        & scftol, iScf, d_spectrum_max, tSpectrumConverged)
+    commutator_max = maxval(abs(commutator))
+    tCommutatorConverged = (maxval(abs(commutator)) < scftol)
 
-    write(*, '(I4,2X,3(1X,F16.9),3X,E16.9)') iScf, total_ene, exchange_energy, x_en_2, change_max
+    write(*, '(I4,2X,3(1X,E16.9),3X,E16.9)') iScf, total_ene, d_spectrum_max, commutator_max, d_pot_max
+
+    ! In case of DIIS, use commutator instead of potential change for convergence
+    if (mixnr == 3) then
+      tConverged = tCommutatorConverged
+    end if
 
     ! if self-consistency is reached, exit loop
     if (tConverged) exit lpScf
@@ -233,7 +245,7 @@ program HFAtom
         & exchange_energy, x_en_2, conf_energy, total_ene)
   end if
 
-  write(*, '(A,E20.12)') 'Potential Matrix Elements converged to ', change_max
+  write(*, '(A,E20.12)') 'Potential Matrix Elements converged to ', d_pot_max
   write(*, '(A)') ' '
 
   if (tZora) then