diff --git a/io/cdn_io.F90 b/io/cdn_io.F90
index 962a82a9036e63f6e137e9a6608eeb8ca3230597..92941fb3ae848e5a46c0d87f385c5ec73ad5b4cc 100644
--- a/io/cdn_io.F90
+++ b/io/cdn_io.F90
@@ -20,6 +20,7 @@ MODULE m_cdn_io
    USE m_wrtdop
    USE m_cdnpot_io_hdf
    USE m_cdnpot_io_common
+   USE m_constants
 #ifdef CPP_HDF
    USE hdf5
 #endif
@@ -29,7 +30,9 @@ MODULE m_cdn_io
    PUBLIC printDensityFileInfo
    PUBLIC readDensity, writeDensity
    PUBLIC isDensityFilePresent, isCoreDensityPresent
+   PUBLIC isDensityMatrixPresent
    PUBLIC readCoreDensity, writeCoreDensity
+   PUBLIC readDensityMatrix
    PUBLIC readStars, writeStars
    PUBLIC readStepfunction, writeStepfunction
    PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
@@ -1154,6 +1157,60 @@ MODULE m_cdn_io
 
    END SUBROUTINE readStepfunction
 
+   SUBROUTINE readDensityMatrix(input, atoms, n_mmp, l_error)
+
+      TYPE(t_input), INTENT(IN)  :: input
+      TYPE(t_atoms), INTENT(IN)  :: atoms
+      COMPLEX, INTENT(OUT)       :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
+                                          -lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
+                                          atoms%n_u,input%jspins)
+      LOGICAL, INTENT(OUT)       :: l_error
+
+      INTEGER                    :: mode, ioStatus
+      LOGICAL                    :: l_exist
+      CHARACTER(LEN=30)          :: filename
+
+      l_error = .FALSE.
+      ioStatus = 0
+      CALL getMode(mode)
+
+      IF(mode.EQ.CDN_HDF5_MODE) THEN
+         l_exist = .FALSE.
+         IF(.NOT.l_exist) THEN
+            mode = CDN_STREAM_MODE
+         END IF
+      END IF
+
+      IF(mode.EQ.CDN_STREAM_MODE) THEN
+         INQUIRE(FILE='cdn.str',EXIST=l_exist)
+         IF (l_exist) THEN
+            STOP 'cdn.str code path not yet implemented!'
+         END IF
+         IF (.NOT.l_exist) THEN
+            mode = CDN_DIRECT_MODE
+         END IF
+      END IF
+
+      IF (mode.EQ.CDN_DIRECT_MODE) THEN
+         filename = 'n_mmp_mat'
+         INQUIRE(FILE=TRIM(ADJUSTL(filename)),EXIST=l_exist)
+         IF(.NOT.l_exist) THEN
+            l_error = .TRUE.
+            RETURN
+         END IF
+
+         OPEN (69,file=TRIM(ADJUSTL(filename)),status='unknown',form='formatted')
+         READ (69,'(7f20.13)',IOSTAT=ioStatus) n_mmp
+         CLOSE(69)
+!         IF (ioStatus.NE.0) THEN
+!            l_error = .TRUE.
+!            RETURN
+!         END IF
+
+      END IF
+   
+   END SUBROUTINE readDensityMatrix
+
    SUBROUTINE setStartingDensity(l_noco)
 
       LOGICAL,INTENT(IN) :: l_noco
@@ -1431,5 +1488,38 @@ MODULE m_cdn_io
       isCoreDensityPresent = .FALSE.
    END FUNCTION isCoreDensityPresent
 
+   LOGICAL FUNCTION isDensityMatrixPresent()
+
+      LOGICAL             :: l_exist
+      INTEGER             :: mode
+      CHARACTER(LEN=12)   :: filename
+
+      CALL getMode(mode)
+
+      filename="n_mmp_mat"
+      INQUIRE (file=filename,exist=l_exist)
+      IF(l_exist) THEN
+         isDensityMatrixPresent = l_exist
+         RETURN
+      END IF
+
+!      IF (mode.EQ.CDN_STREAM_MODE) THEN
+!         STOP 'Not yet implemented! (CDN_STREAM_MODE in isDensityMatrixPresent)'
+!      END IF
+
+      IF (mode.EQ.CDN_HDF5_MODE) THEN
+!#ifdef CPP_HDF
+!         STOP 'Not yet implemented (uncomment line)!'
+!         l_exist = isDensityMatrixPresentHDF()
+!         IF(l_exist) THEN
+!            isDensityMatrixPresent = l_exist
+!            RETURN
+!         END IF
+!#endif
+      END IF
+
+      isDensityMatrixPresent = .FALSE.
+
+   END FUNCTION isDensityMatrixPresent
 
 END MODULE m_cdn_io
diff --git a/ldau/u_setup.f90 b/ldau/u_setup.f90
index 8649252177d139f312c428c10f05cafd9bc48ce8..b03fe5acd376705826f2d3e45a697fb28d48738f 100644
--- a/ldau/u_setup.f90
+++ b/ldau/u_setup.f90
@@ -24,6 +24,7 @@ CONTAINS
     USE m_vmmp
     USE m_types
     USE m_constants
+    USE m_cdn_io
     IMPLICIT NONE
     TYPE(t_sym),INTENT(IN)          :: sym
     TYPE(t_results),INTENT(INOUT)   :: results
@@ -38,7 +39,7 @@ CONTAINS
     INTEGER itype,ispin,j,k,l,jspin,urec,i_u
     INTEGER noded,nodeu,ios,lty(atoms%n_u)
     REAL wronk
-    LOGICAL n_mmp_exist,n_exist
+    LOGICAL n_mmp_exist,n_exist, l_error
     CHARACTER*8 l_type*2,l_form*9
     CHARACTER*12 ::filename
     REAL f(atoms%jmtd,2),g(atoms%jmtd,2),theta(atoms%n_u),phi(atoms%n_u),zero(atoms%n_u)
@@ -48,9 +49,7 @@ CONTAINS
     !
     ! look, whether density matrix exists already:
     !
-    filename="n_mmp_mat"
-    INQUIRE (file=filename,exist=n_mmp_exist)
-    IF (n_mmp_exist.AND.atoms%n_u>0) THEN
+    IF (isDensityMatrixPresent().AND.atoms%n_u>0) THEN
        !
        ! calculate slater integrals from u and j
        !
@@ -74,10 +73,8 @@ CONTAINS
        ! read density matrix
        !
        ALLOCATE (ns_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u,input%jspins))
-       OPEN (69,file=filename,status='unknown',form='formatted')
-       READ (69,9000) ns_mmp
-9000   FORMAT(7f20.13)
-       CLOSE (69)
+       CALL readDensityMatrix(input,atoms,ns_mmp,l_error)
+       IF(l_error) CALL juDFT_error('Error in reading density matrix!',calledby='u_setup')
        !
        ! check for possible rotation of n_mmp
        !