diff --git a/common/lib/accuracy.F90 b/common/lib/accuracy.F90
index 19ce9a94..2d2fc8bf 100644
--- a/common/lib/accuracy.F90
+++ b/common/lib/accuracy.F90
@@ -3,7 +3,7 @@
 !! Not all routines use the string length specifications to set their character string lengths.
 module common_accuracy
 
-  use, intrinsic :: iso_fortran_env, only : real64
+  use, intrinsic :: iso_fortran_env, only : real64, real128
 
   implicit none
   private
diff --git a/slateratom/lib/diagonalizations.f90 b/slateratom/lib/diagonalizations.f90
index b488c319..2703ae7f 100644
--- a/slateratom/lib/diagonalizations.f90
+++ b/slateratom/lib/diagonalizations.f90
@@ -34,8 +34,9 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
     !! eigenvalues of overlap matrices
     real(dp), allocatable :: eigenvalues(:)
 
-    !> auxiliary variables
+    !! auxiliary variables
     integer :: ll, diagsize, ii
+    real(dp) :: e_min, e_max, cond_num
 
     do ll = 0, max_l
 
@@ -57,11 +58,30 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
         stop
       end if
 
+      if (ll == 0) then
+        e_min = minval(eigenvalues)
+        e_max = maxval(eigenvalues)
+      end if
+
+      if (minval(eigenvalues) < e_min) then
+        e_min = minval(eigenvalues)
+      end if
+
+      if (maxval(eigenvalues) > e_max) then
+        e_max = maxval(eigenvalues)
+      end if
+
       deallocate(overlap, eigenvalues)
 
     end do
+
+    cond_num = abs(e_max / e_min)
+    write(*, '(A,E16.8)') 'Condition number is ', cond_num
+    ! 16 is due to double precision
+    write(*, '(A,I2)') 'Expect convergence only up to approx. 1e-', 16 - int(anint(log10(cond_num)))
     write(*,*) ' '
 
+
   end subroutine diagonalize_overlap
 
 
diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90
index 042fd573..3b53cc09 100644
--- a/slateratom/lib/globals.f90
+++ b/slateratom/lib/globals.f90
@@ -85,7 +85,10 @@ module globals
   real(dp) :: commutator_max
 
   !> density matrix supervector
-  real(dp), allocatable :: pp(:,:,:,:)
+  real(dp), allocatable :: pp(:,:,:,:), pp_old(:,:,:,:)
+
+  !> density matrix max difference
+  real(dp) :: pp_diff
 
   !> fock matrix supervector
   real(dp), allocatable :: ff(:,:,:,:)
@@ -180,8 +183,8 @@ module globals
   !> true, if zero-order regular approximation for relativistic effects is desired
   logical :: tZora
 
-  !> true, if SCF cycle reached convergency on a quantity
-  logical :: tCommutatorConverged, tEnergyConverged, tEigenspectrumConverged
+  !> true, if SCF cycle reached convergency on a given quantity
+  logical :: tCommutatorConverged, tEnergyConverged, tEigenspectrumConverged, tDensityConverged
 
   !> identifier of mixer
   integer :: mixnr
@@ -262,6 +265,7 @@ subroutine allocate_globals()
     ! fourth index of cof is the eigenvalue index, see densmatrix build
     allocate(cof(2, 0:max_l, problemsize, problemsize))
     allocate(pp(2, 0:max_l, problemsize, problemsize))
+    allocate(pp_old(2, 0:max_l, problemsize, problemsize))
 
     weight(:) = 0.0_dp
     abcissa(:) = 0.0_dp
diff --git a/slateratom/lib/utilities.f90 b/slateratom/lib/utilities.f90
index f01a26f1..08536702 100644
--- a/slateratom/lib/utilities.f90
+++ b/slateratom/lib/utilities.f90
@@ -7,7 +7,8 @@ module utilities
   private
 
   public :: check_electron_number, check_convergence_energy, check_convergence_commutator,&
-      & check_convergence_eigenspectrum, compute_commutator, check_convergence_pot
+      & check_convergence_eigenspectrum, compute_commutator, check_convergence_pot,&
+      & check_convergence_density
   public :: vector_length, fak, zeroOutCpotOfEmptyDensitySpinChannels
 
 
@@ -164,7 +165,37 @@ subroutine compute_commutator(max_l, num_alpha, poly_order, ff, pp, ss, commutat
   end subroutine 
 
 
-  !> Checks convergence by computing energy change from the last SCF iteration.
+  !> Checks convergence by computing max. density change from the last SCF iteration.
+  pure subroutine check_convergence_density(pp_old, pp_new, scftol, iScf, res, tConverged)
+
+    !> old and new energy to compare
+    real(dp), intent(in) :: pp_old(:,:,:,:), pp_new(:,:,:,:)
+
+    !> scf tolerance, i.e. convergence criteria
+    real(dp), intent(in) :: scftol
+
+    !> current SCF step
+    integer, intent(in) :: iScf
+
+    !> obtained maximum change in energy
+    real(dp), intent(out) :: res
+
+    !> true, if SCF converged
+    logical, intent(out) :: tConverged
+
+    res = 0.0_dp
+
+    res = maxval(abs(pp_new - pp_old))
+    tConverged = res < scftol
+
+    if (iScf < 3) then
+      tConverged = .false.
+    end if
+
+  end subroutine check_convergence_density
+
+
+    !> Checks convergence by computing energy change from the last SCF iteration.
   pure subroutine check_convergence_energy(energy_old, energy_new, scftol, iScf, change,&
       & tConverged)
 
@@ -197,7 +228,7 @@ end subroutine check_convergence_energy
 
   !> Checks convergence by evaluating change in the occupied part of the eigenspectrum
   pure subroutine check_convergence_eigenspectrum(max_l, num_alpha, poly_order, eigval_new,&
-      & eigval_old, occ, scftol, iScf, change, tConverged)
+      & eigval_old, occ, scftol, iScf, res, tConverged)
 
     !> maximum angular momentum
     integer, intent(in) :: max_l
@@ -221,7 +252,7 @@ pure subroutine check_convergence_eigenspectrum(max_l, num_alpha, poly_order, ei
     integer, intent(in) :: iScf
 
     !> obtained change
-    real(dp), intent(out) :: change
+    real(dp), intent(out) :: res
 
     !> true, if SCF converged
     logical, intent(out) :: tConverged
@@ -232,22 +263,28 @@ pure subroutine check_convergence_eigenspectrum(max_l, num_alpha, poly_order, ei
     ! Occupation
     real(dp) :: occ_ii
 
-    change = 0.0_dp
+    ! Difference
+    real(dp) :: e_diff
 
-    ! Frobenius norm, but only occupied states contribute 
+    ! 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
           occ_ii = abs(occ(iSpin, ll, ii))
           if (occ_ii > 1e-16) then
-            change = change + (eigval_new(iSpin, ll, ii) - eigval_old(iSpin, ll, ii))**2 
+            ! change = change + (eigval_new(iSpin, ll, ii) - eigval_old(iSpin, ll, ii))**2 
+            e_diff = abs(eigval_new(iSpin, ll, ii) - eigval_old(iSpin, ll, ii)) 
+            if (e_diff > res) then
+              res = e_diff
+            end if
           end if
         end do
       end do
     end do
 
-    tConverged = change < scftol
+    tConverged = res < scftol
 
     if (iScf < 3) then
       tConverged = .false.
diff --git a/slateratom/prog/main.f90 b/slateratom/prog/main.f90
index 8b476cee..9448c686 100644
--- a/slateratom/prog/main.f90
+++ b/slateratom/prog/main.f90
@@ -17,7 +17,7 @@ program HFAtom
   use dft, only : check_accuracy, density_grid
   use utilities, only : check_electron_number, check_convergence_energy,&
       & check_convergence_commutator, check_convergence_eigenspectrum,&
-      & compute_commutator, check_convergence_pot
+      & compute_commutator, check_convergence_density
   use zora_routines, only : scaled_zora
   use cmdargs, only : parse_command_arguments
   use common_poisson, only : TBeckeGridParams
@@ -82,8 +82,7 @@ program HFAtom
   ! WARNING: too high number of grid points somehow
   ! manages to break the Broyden mixer!
   if (xcFunctional%isMGGA(xcnr)) then
-    ! num_mesh_points = num_mesh_points + 2000
-    num_mesh_points = num_mesh_points + 500
+    num_mesh_points = num_mesh_points + 1000
   end if
 
   call echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha, conf_r0,&
@@ -146,18 +145,20 @@ program HFAtom
       & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, vtau, alpha, pot_old,&
       & pot_new, tZora, ff, commutator, camAlpha, camBeta)
 
+  pp_old(:,:,:,:) = pp
+
   ! self-consistency cycles
   write(*,*) 'Energies in Hartree'
   write(*,*)
-  write(*,*) ' Iter |   Total energy  |   HF-X energy  |   XC energy   |   max(abs([F,PS]))   &
-      & | Delta (Total energy) | Delta (spectrum)'
+  write(*,*) ' Iter |   Total energy    |        P          |       [F,PS]      |         E         |       lambda'
   write(*,*) '-------------------------------------------------------------------------------------&
-      & ------------------------------------------'
+      &--------------------'
   lpScf: do iScf = 1, maxiter
 
     pot_old(:,:,:,:) = pot_new
     total_ene_old = total_ene
     eigval_old(:,:,:) = eigval
+    pp_old(:,:,:,:) = pp
 
     ! diagonalize
     call diagonalize(max_l, num_alpha, poly_order, ff, ss, cof, eigval)
@@ -190,16 +191,14 @@ program HFAtom
 
     call check_convergence_commutator(commutator, scftol, iScf, commutator_max,&
         & tCommutatorConverged)
-    ! For debugging:
-    ! call check_convergence_pot(pot_old, pot_new, max_l, problemsize, scftol, iScf,&
-        ! & commutator_max, tCommutatorConverged)
+    call check_convergence_density(pp_old, pp, scftol, iScf, pp_diff, tDensityConverged)
     call check_convergence_energy(total_ene_old, total_ene, scftol, iScf, total_ene_diff,&
         & tEnergyConverged)
     call check_convergence_eigenspectrum(max_l, num_alpha, poly_order, eigval, eigval_old, occ,&
         & scftol, iScf, eigval_diff, tEigenspectrumConverged)
 
     ! Print SCF loop information
-    write(*, '(I4,2X,3(1X,F16.9),7X,E16.9,8X,E16.9,8X,E16.9)') iScf, total_ene, exchange_energy, x_en_2,&
+    write(*, '(I4,5X,E16.9,4X,E16.9,4X,E16.9,4X,E16.9,4X,E16.9)') iScf, total_ene, pp_diff,&
         & commutator_max, total_ene_diff, eigval_diff
 
     ! if self-consistency is reached, exit loop
@@ -246,6 +245,7 @@ program HFAtom
   end if
 
   write(*, '(A,E20.12)') '[F,PS] converged to', commutator_max
+  write(*, '(A,E20.12)') 'Density matrix converged to', pp_diff
   write(*, '(A)') ' '
 
   if (tZora) then