diff --git a/devel/density_matrix/.gitignore b/devel/density_matrix/.gitignore
deleted file mode 100644
index 7ac9fbf..0000000
--- a/devel/density_matrix/.gitignore
+++ /dev/null
@@ -1,5 +0,0 @@
-IRPF90_temp/
-IRPF90_man/
-irpf90.make
-irpf90_entities
-tags
\ No newline at end of file
diff --git a/devel/density_matrix/NEED b/devel/density_matrix/NEED
deleted file mode 100644
index d3d4d2c..0000000
--- a/devel/density_matrix/NEED
+++ /dev/null
@@ -1 +0,0 @@
-determinants
diff --git a/devel/density_matrix/README.rst b/devel/density_matrix/README.rst
deleted file mode 100644
index ded0613..0000000
--- a/devel/density_matrix/README.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-====================
-DensityMatrix Module
-====================
-
-Prints the 1- and 2- body density matrices
diff --git a/devel/density_matrix/density_matrix_array.irp.f b/devel/density_matrix/density_matrix_array.irp.f
deleted file mode 100644
index 51ab45c..0000000
--- a/devel/density_matrix/density_matrix_array.irp.f
+++ /dev/null
@@ -1,115 +0,0 @@
- use bitmasks
- BEGIN_PROVIDER [ double precision, two_e_dm_aa, (mo_num,mo_num,mo_num,mo_num) ]
-&BEGIN_PROVIDER [ double precision, two_e_dm_bb, (mo_num,mo_num,mo_num,mo_num) ]
-&BEGIN_PROVIDER [ double precision, two_e_dm_ab, (mo_num,mo_num,mo_num,mo_num) ]
-   implicit none
-   use bitmasks
-   BEGIN_DOC
-   ! Temporary files for 2-e dm calculation
-   END_DOC
-   integer                        :: getUnitAndOpen
-   
-   integer                        :: k,l,degree, idx,i,j
-   integer                        :: exc(0:2,2,2),n_occ_alpha
-   double precision               :: phase, coef
-   integer                        :: h1,h2,p1,p2,s1,s2, e1, e2
-   double precision               :: ck, cl
-   character*(128), parameter     :: f = '(i8,4(x,i5),x,d16.8)'
-   integer :: istate
-
-   two_e_dm_aa = 0.d0
-   two_e_dm_ab = 0.d0
-   two_e_dm_bb = 0.d0
-
-   istate = 1
-   ! OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,ck,ckl,i,j,e1,e2,cl,phase,h1,p1,h2,p2,s1,s2,occ)
-   ! OMP DO SCHEDULE(dynamic,64)
-   do k=1,N_det
-     ck = psi_coef(k,istate)
-     call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
-     call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
-     ckl = psi_coef(k,istate) * psi_coef(k,istate)
-     do i = 1,elec_alpha_num
-       e1=occ(i,1)
-       do j = 1,elec_alpha_num
-         e2=occ(j,1)
-         ! alpha-alpha
-         two_e_dm_aa(e1,e2,e1,e2) += 0.5d0*ckl
-         two_e_dm_aa(e1,e2,e2,e1) -= 0.5d0*ckl
-       enddo
-       do j = 1,elec_beta_num
-         e2=occ(j,2)
-         ! alpha-beta
-         two_e_dm_ab(e1,e2,e1,e2) += ckl
-       enddo
-     enddo
-     do i = 1,elec_beta_num
-       e1=occ(i,2)
-       do j = 1,elec_beta_num
-         e2=occ(j,2)
-         ! beta-beta
-         two_e_dm_bb(e1,e2,e1,e2) += 0.5d0*ckl
-         two_e_dm_bb(e1,e2,e2,e1) -= 0.5d0*ckl
-       enddo
-     enddo
-     
-     do l=1,k-1
-       cl = 2.d0*psi_coef(l,istate)
-       call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
-       if (degree == 2) then
-         call get_double_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
-         call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-         ckl = phase*ck*cl
-         select case (s1+s2)
-           case(2) ! alpha alpha
-             two_e_dm_aa(h1,h2,p1,p2) += ckl
-             two_e_dm_aa(h1,h2,p2,p1) -= ckl
-           case(3) ! alpha beta
-             two_e_dm_ab(h1,h2,p1,p2) += ckl
-           case(4) ! beta beta
-             two_e_dm_bb(h1,h2,p1,p2) += ckl
-             two_e_dm_bb(h1,h2,p2,p1) -= ckl
-         end select
-       else if (degree == 1) then
-         call get_single_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
-         call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-         double precision               :: ckl
-         ckl = phase*ck*cl
-         call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
-         call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
-         select case (s1)
-           case (1)  ! Alpha single excitation
-             integer                        :: occ(N_int*bit_kind_size,2)
-             do i  = 1, elec_alpha_num
-               p2=occ(i,1)
-               h2=p2
-               two_e_dm_aa(h1,h2,p1,p2) += ckl
-               two_e_dm_aa(h1,h2,p2,p1) -= ckl
-             enddo
-             do i = 1, elec_beta_num
-               p2=occ(i,2)
-               h2=p2
-               two_e_dm_ab(h1,h2,p1,p2) += ckl
-             enddo
-           case (2)  ! Beta single excitation
-             do i  = 1, elec_alpha_num
-               p2=occ(i,1)
-               h2=p2
-               two_e_dm_ab(h1,h2,p1,p2) += ckl
-             enddo
-             do i = 1, elec_beta_num
-               p2=occ(i,2)
-               h2=p2
-               two_e_dm_bb(h1,h2,p1,p2) += ckl
-               two_e_dm_bb(h1,h2,p2,p1) -= ckl
-             enddo
-         end select
-       endif
-     enddo
-   enddo
-   ! OMP END DO
-   ! OMP END PARALLEL
-   
-END_PROVIDER
- 
- 
diff --git a/devel/density_matrix/print_2rdm.irp.f b/devel/density_matrix/print_2rdm.irp.f
deleted file mode 100644
index 682222a..0000000
--- a/devel/density_matrix/print_2rdm.irp.f
+++ /dev/null
@@ -1,43 +0,0 @@
-program print_2rdm
-  implicit none
-  integer                        :: i,j,k,l
-  double precision, external     :: get_two_e_integral
-  read_wf = .True.
-  TOUCH read_wf
-  
-  double precision               :: e(10)
-  double precision, parameter    :: thr = 1.d-15
-  e = 0.d0
-  
-  print *,  '1RDM'
-  do i=1,mo_num
-    do j=1,mo_num
-      if (dabs(one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)) > thr) then
-        print *,  i, j, one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)
-      endif
-      e(4) += one_e_dm_mo_alpha(i,j,1) * mo_one_e_integrals(i,j)
-      e(4) += one_e_dm_mo_beta(i,j,1) * mo_one_e_integrals(i,j)
-    enddo
-  enddo
-  
-  print *,  '2RDM'
-  do i=1,mo_num
-    do j=1,mo_num
-      do k=1,mo_num
-        do l=1,mo_num
-          if (dabs(two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)) > thr) then
-            print *,  i, j, k, l, two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)
-          endif
-          e(1) += two_e_dm_aa(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-          e(2) += two_e_dm_bb(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-          e(3) += two_e_dm_ab(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-        enddo
-      enddo
-    enddo
-  enddo
-  
-  print *,  ''
-  print *,  'Energy ', sum(e(1:4)) + nuclear_repulsion
-  
-  
-end
diff --git a/devel/density_matrix/print_2rdm_decomposed.irp.f b/devel/density_matrix/print_2rdm_decomposed.irp.f
deleted file mode 100644
index 1b94fe5..0000000
--- a/devel/density_matrix/print_2rdm_decomposed.irp.f
+++ /dev/null
@@ -1,78 +0,0 @@
-program print_2rdm_decomposed
-  implicit none
-  integer                        :: i,j,k,l
-  double precision, external     :: get_two_e_integral
-  read_wf = .True.
-  TOUCH read_wf
-  
-  double precision               :: e(10)
-  double precision, parameter    :: thr = 1.d-15
-  e = 0.d0
-  
-  print *,  '1RDM ALPHA'
-  do i=1,mo_num
-    do j=1,mo_num
-      if (dabs(one_e_dm_mo_alpha(i,j,1)) > thr) then
-        print *,  i, j, one_e_dm_mo_alpha(i,j,1)
-      endif
-      e(4) += one_e_dm_mo_alpha(i,j,1) * mo_one_e_integrals(i,j)
-    enddo
-  enddo
-  
-  print *,  '1RDM BETA'
-  do i=1,mo_num
-    do j=1,mo_num
-      if (dabs(one_e_dm_mo_beta(i,j,1)) > thr) then
-        print *,  i, j, one_e_dm_mo_beta(i,j,1)
-      endif
-      e(4) += one_e_dm_mo_beta(i,j,1) * mo_one_e_integrals(i,j)
-    enddo
-  enddo
-
- print *,  '2RDM ALPHA ALPHA'
- do i=1,mo_num
-   do j=1,mo_num
-     do k=1,mo_num
-       do l=1,mo_num
-         if (dabs(two_e_dm_aa(i,j,k,l)) > thr) then
-          print *,  i, j, k, l, two_e_dm_aa(i,j,k,l)
-         endif
-         e(1) += two_e_dm_aa(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-       enddo
-     enddo
-   enddo
- enddo
-
- print *,  '2RDM BETA BETA'
- do i=1,mo_num
-   do j=1,mo_num
-     do k=1,mo_num
-       do l=1,mo_num
-         if (dabs(two_e_dm_bb(i,j,k,l)) > thr) then
-          print *,  i, j, k, l, two_e_dm_bb(i,j,k,l)
-         endif
-         e(2) += two_e_dm_bb(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-       enddo
-     enddo
-   enddo
- enddo
-
- print *,  '2RDM ALPHA BETA'
- do i=1,mo_num
-   do j=1,mo_num
-     do k=1,mo_num
-       do l=1,mo_num
-         if (dabs(two_e_dm_ab(i,j,k,l)) > thr) then
-          print *,  i, j, k, l, two_e_dm_ab(i,j,k,l)
-         endif
-         e(3) += two_e_dm_ab(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
-       enddo
-     enddo
-   enddo
- enddo
-
- print *,  ''
- print *,  'Energy ', sum(e(1:4)) + nuclear_repulsion
-
-
-end
diff --git a/devel/import_integrals_ao/.gitignore b/devel/import_integrals_ao/.gitignore
new file mode 100644
index 0000000..1561915
--- /dev/null
+++ b/devel/import_integrals_ao/.gitignore
@@ -0,0 +1,59 @@
+IRPF90_temp/
+IRPF90_man/
+build.ninja
+irpf90.make
+ezfio_interface.irp.f
+irpf90_entities
+tags
+Makefile
+ao_basis
+ao_one_e_ints
+ao_two_e_erf_ints
+ao_two_e_ints
+aux_quantities
+becke_numerical_grid
+bitmask
+cis
+cisd
+cipsi
+davidson
+davidson_dressed
+davidson_undressed
+density_for_dft
+determinants
+dft_keywords
+dft_utils_in_r
+dft_utils_one_e
+dft_utils_two_body
+dressing
+dummy
+electrons
+ezfio_files
+fci
+generators_cas
+generators_full
+hartree_fock
+iterations
+kohn_sham
+kohn_sham_rs
+mo_basis
+mo_guess
+mo_one_e_ints
+mo_two_e_erf_ints
+mo_two_e_ints
+mpi
+mrpt_utils
+nuclei
+perturbation
+pseudo
+psiref_cas
+psiref_utils
+scf_utils
+selectors_cassd
+selectors_full
+selectors_utils
+single_ref_method
+slave
+tools
+utils
+zmq
diff --git a/devel/import_integrals_ao/NEED b/devel/import_integrals_ao/NEED
new file mode 100644
index 0000000..2637b96
--- /dev/null
+++ b/devel/import_integrals_ao/NEED
@@ -0,0 +1 @@
+nuclei ao_one_e_ints ao_two_e_ints
diff --git a/devel/import_integrals_ao/README.rst b/devel/import_integrals_ao/README.rst
new file mode 100644
index 0000000..e96347c
--- /dev/null
+++ b/devel/import_integrals_ao/README.rst
@@ -0,0 +1,28 @@
+===================
+import_integrals_ao
+===================
+
+Module to read all the integrals in the |AO| basis from files (all in atomic units).
+
+
+The following files are required:
+
+- :file:`S.qp` : overlap integrals
+- :file:`T.qp` : kinetic integrals
+- :file:`V.qp` : electron-nucleus potential integrals
+- :file:`P.qp` : pseudo-potential integrals
+- :file:`W.qp` : electron repulsion integrals
+
+If present, the :file:`E.qp` file, should contain the nuclear repulsion energy.
+
+In all the other files, there is one integral per line and for the one-electron integral
+$\int \chi_i(r) \hat{O} \chi_j(r) dr$, the format is
+
+  i j value
+
+and for two electron integral the format uses the physicists' convention,
+$\int \chi_i(r_1) \chi_j(r_2) \hat{O} \chi_k(r_1) \chi_l(r_2) dr$:
+
+  i j k l value
+
+
diff --git a/devel/import_integrals_ao/import_integrals_ao.irp.f b/devel/import_integrals_ao/import_integrals_ao.irp.f
new file mode 100644
index 0000000..fb7a638
--- /dev/null
+++ b/devel/import_integrals_ao/import_integrals_ao.irp.f
@@ -0,0 +1,108 @@
+program print_integrals
+  print *, 'Number of AOs?'
+  read(*,*) ao_num
+  TOUCH ao_num
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j,k,l
+  double precision :: integral
+  double precision, allocatable :: A(:,:)
+
+  integer             :: n_integrals 
+  integer(key_kind), allocatable   :: buffer_i(:) 
+  real(integral_kind), allocatable :: buffer_values(:)
+
+  call ezfio_set_ao_basis_ao_num(ao_num)
+
+  allocate (A(ao_num,ao_num))
+  
+  A(1,1) = huge(1.d0)
+  iunit = getunitandopen('E.qp','r')
+  read (iunit,*,end=9) A(1,1)
+  9 continue
+  close(iunit)
+  if (A(1,1) /= huge(1.d0)) then
+    call ezfio_set_nuclei_nuclear_repulsion(A(1,1))
+    call ezfio_set_nuclei_io_nuclear_repulsion("Read")
+  endif
+
+  A = 0.d0
+  iunit = getunitandopen('T.qp','r')
+  do 
+    read (iunit,*,end=10) i,j, integral
+    A(i,j) = integral
+  enddo
+  10 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read")
+
+  A = 0.d0
+  iunit = getunitandopen('S.qp','r')
+  do 
+    read (iunit,*,end=11) i,j, integral
+    A(i,j) = integral
+  enddo
+  11 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read")
+
+  A = 0.d0
+  iunit = getunitandopen('P.qp','r')
+  do 
+    read (iunit,*,end=14) i,j, integral
+    A(i,j) = integral
+  enddo
+  14 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read")
+
+  A = 0.d0
+  iunit = getunitandopen('V.qp','r')
+  do 
+    read (iunit,*,end=12) i,j, integral
+    A(i,j) = integral
+  enddo
+  12 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_e_n(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_e_n("Read")
+
+  allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
+  iunit = getunitandopen('W.qp','r')
+  n_integrals=0
+  buffer_values = 0.d0
+  do 
+    read (iunit,*,end=13) i,j,k,l, integral
+    n_integrals += 1
+    call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) )
+    buffer_values(n_integrals) = integral
+    if (n_integrals == size(buffer_i)) then
+      call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
+      n_integrals = 0
+    endif
+  enddo
+  13 continue
+  close(iunit)
+  
+  if (n_integrals > 0) then
+    call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
+  endif
+
+  call map_sort(ao_integrals_map)
+  call map_unique(ao_integrals_map)
+
+  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+  call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+  
+end