diff --git a/plugins/DensityMatrix/NEEDED_CHILDREN_MODULES b/plugins/DensityMatrix/NEEDED_CHILDREN_MODULES
index a953d6a6..065099eb 100644
--- a/plugins/DensityMatrix/NEEDED_CHILDREN_MODULES
+++ b/plugins/DensityMatrix/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Bitmask
\ No newline at end of file
+Bitmask Determinants
diff --git a/plugins/DensityMatrix/density_matrix.irp.f b/plugins/DensityMatrix/density_matrix.irp.f
deleted file mode 100644
index 5e7d7cec..00000000
--- a/plugins/DensityMatrix/density_matrix.irp.f
+++ /dev/null
@@ -1,211 +0,0 @@
- use bitmasks
- BEGIN_PROVIDER [ integer, iunit_two_body_dm_aa ]
-&BEGIN_PROVIDER [ integer, iunit_two_body_dm_ab ]
-&BEGIN_PROVIDER [ integer, iunit_two_body_dm_bb ]
-   implicit none
-   use bitmasks
-   BEGIN_DOC
-   ! Temporary files for 2-body dm calculation
-   END_DOC
-   integer                        :: getUnitAndOpen
-   
-   iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','w')
-   iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','w')
-   iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','w')
-   ! Compute two body DM in file
-   integer                        :: k,l,degree, idx,i
-   integer                        :: exc(0:2,2,2),n_occ_alpha
-   double precision               :: phase, coef
-   integer                        :: h1,h2,p1,p2,s1,s2
-   double precision               :: ck, cl
-   character*(128), parameter     :: f = '(i8,4(x,i5),x,d16.8)'
-   do k=1,det_num
-     ck = (det_coef_provider(k)+det_coef_provider(k))
-     do l=1,k-1
-       cl = det_coef_provider(l)
-       call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int)
-       if (degree == 2) then
-         call get_double_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
-         call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-         call bielec_integrals_index(h1,h2,p1,p2,idx)
-         ckl = phase*ck*cl
-         select case (s1+s2)
-           case(2) ! alpha alpha
-             write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
-             call bielec_integrals_index(h1,h2,p2,p1,idx)
-             write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
-           case(3) ! alpha beta
-             write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
-           case(4) ! beta beta
-             write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
-             call bielec_integrals_index(h1,h2,p2,p1,idx)
-             write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
-         end select
-       else if (degree == 1) then
-         call get_mono_excitation(det_provider(1,1,k),det_provider(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(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
-         call bitstring_to_list(det_provider(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
-               call bielec_integrals_index(h1,h2,p1,p2,idx)
-               write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
-               call bielec_integrals_index(h1,h2,p2,p1,idx)
-               write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
-             enddo
-             do i = 1, elec_beta_num
-               p2=occ(i,2)
-               h2=p2
-               call bielec_integrals_index(h1,h2,p1,p2,idx)
-               write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
-             enddo
-           case (2)  ! Beta single excitation
-             do i  = 1, elec_alpha_num
-               p2=occ(i,1)
-               h2=p2
-               call bielec_integrals_index(h1,h2,p1,p2,idx)
-               write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
-             enddo
-             do i = 1, elec_beta_num
-               p2=occ(i,2)
-               h2=p2
-               call bielec_integrals_index(h1,h2,p1,p2,idx)
-               write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
-               call bielec_integrals_index(h1,h2,p2,p1,idx)
-               write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
-             enddo
-         end select
-       endif
-     enddo
-   enddo
-   ! Sort file
-   ! Merge coefs
-   
-   close(iunit_two_body_dm_aa)
-   close(iunit_two_body_dm_ab)
-   close(iunit_two_body_dm_bb)
-   character*(128)                :: filename
-   filename = trim(ezfio_filename)//'/work/two_body_aa.tmp'
-   call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
-   filename = trim(ezfio_filename)//'/work/two_body_ab.tmp'
-   call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
-   filename = trim(ezfio_filename)//'/work/two_body_bb.tmp'
-   call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
-   iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','r')
-   iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','r')
-   iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','r')
-END_PROVIDER
- 
- 
- BEGIN_TEMPLATE
- 
-BEGIN_PROVIDER [ integer, size_two_body_dm_$AA ]
-   implicit none
-   use bitmasks
-   BEGIN_DOC
-   ! Size of the two body $ALPHA density matrix
-   END_DOC
-   integer *8                     :: key, key_old
-   rewind(iunit_two_body_dm_$AA)
-   size_two_body_dm_$AA = 0
-   key     = 0_8
-   key_old = key
-   do while (.True.)
-     read(iunit_two_body_dm_$AA,*,END=99) key
-     if (key /= key_old) then
-       size_two_body_dm_$AA += 1
-       key_old = key
-     endif
-   end do
-   99 continue
-   
-END_PROVIDER
- 
- BEGIN_PROVIDER [ integer, two_body_dm_index_$AA, (4,size_two_body_dm_$AA) ]
-&BEGIN_PROVIDER [ double precision, two_body_dm_value_$AA, (size_two_body_dm_$AA) ]
-   implicit none
-   use bitmasks
-   BEGIN_DOC
-   ! Two body $ALPHA density matrix
-   END_DOC
-   rewind(iunit_two_body_dm_$AA)
-   integer *8                     :: key, key_old
-   integer                        :: ii, i,j,k,l
-   double precision               :: c
-   key     = 0_8
-   key_old = key
-   ii = 0
-   do while (.True.)
-     read(iunit_two_body_dm_$AA,*,END=99) key, i,j,k,l, c
-     if (key /= key_old) then
-       ii += 1
-       two_body_dm_index_$AA(1,ii) = i
-       two_body_dm_index_$AA(2,ii) = j
-       two_body_dm_index_$AA(3,ii) = k
-       two_body_dm_index_$AA(4,ii) = l
-       two_body_dm_value_$AA(ii) = 0.d0
-       key_old = key
-     endif
-     two_body_dm_value_$AA(ii) += c
-   enddo
-   99 continue
-   close(iunit_two_body_dm_$AA, status='DELETE')
-END_PROVIDER
- 
- SUBST [ AA, ALPHA ]
- 
- aa ; alpha-alpha ;;
- ab ; alpha-beta  ;;
- bb ; beta-beta   ;;
- 
- END_TEMPLATE
- 
- 
- BEGIN_PROVIDER [ double precision, two_body_dm_diag_aa, (mo_tot_num,mo_tot_num)]
-&BEGIN_PROVIDER [ double precision, two_body_dm_diag_bb, (mo_tot_num,mo_tot_num)]
-&BEGIN_PROVIDER [ double precision, two_body_dm_diag_ab, (mo_tot_num,mo_tot_num)]
-   implicit none
-   use bitmasks
-   BEGIN_DOC
-   ! diagonal part of the two body density matrix
-   END_DOC
-   integer                        :: i,j,k,e1,e2
-   integer                        :: occ(N_int*bit_kind_size,2)
-   double precision               :: ck
-   integer                        :: n_occ_alpha
-   two_body_dm_diag_aa=0.d0
-   two_body_dm_diag_ab=0.d0
-   two_body_dm_diag_bb=0.d0
-   do k = 1, det_num
-     call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
-     call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
-     ck = det_coef_provider(k) * det_coef_provider(k)
-     do i = 1,elec_alpha_num
-       e1=occ(i,1)
-       do j = 1,elec_alpha_num
-         e2=occ(j,1)
-         ! alpha-alpha
-         two_body_dm_diag_aa(e1,e2) = two_body_dm_diag_aa(e1,e2) + ck
-       enddo
-       do j = 1,elec_beta_num
-         e2=occ(j,2)
-         ! alpha-beta
-         two_body_dm_diag_ab(e1,e2) = two_body_dm_diag_ab(e1,e2) + ck
-       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_body_dm_diag_bb(e1,e2) = two_body_dm_diag_bb(e1,e2) + ck
-       enddo
-     enddo
-   enddo
-END_PROVIDER
diff --git a/plugins/DensityMatrix/density_matrix_array.irp.f b/plugins/DensityMatrix/density_matrix_array.irp.f
new file mode 100644
index 00000000..902eac7b
--- /dev/null
+++ b/plugins/DensityMatrix/density_matrix_array.irp.f
@@ -0,0 +1,116 @@
+ use bitmasks
+ BEGIN_PROVIDER [ double precision, two_body_dm_aa, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, two_body_dm_bb, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, two_body_dm_ab, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
+   implicit none
+   use bitmasks
+   BEGIN_DOC
+   ! Temporary files for 2-body dm calculation
+   END_DOC
+   integer                        :: getUnitAndOpen
+   
+   ! Compute two body DM in file
+   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_body_dm_aa = 0.d0
+   two_body_dm_ab = 0.d0
+   two_body_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_body_dm_aa(e1,e2,e1,e2) += 0.5d0*ckl
+         two_body_dm_aa(e1,e2,e2,e1) -= 0.5d0*ckl
+       enddo
+       do j = 1,elec_beta_num
+         e2=occ(j,2)
+         ! alpha-beta
+         two_body_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_body_dm_bb(e1,e2,e1,e2) += 0.5d0*ckl
+         two_body_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_body_dm_aa(h1,h2,p1,p2) += ckl
+             two_body_dm_aa(h1,h2,p2,p1) -= ckl
+           case(3) ! alpha beta
+             two_body_dm_ab(h1,h2,p1,p2) += ckl
+           case(4) ! beta beta
+             two_body_dm_bb(h1,h2,p1,p2) += ckl
+             two_body_dm_bb(h1,h2,p2,p1) -= ckl
+         end select
+       else if (degree == 1) then
+         call get_mono_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_body_dm_aa(h1,h2,p1,p2) += ckl
+               two_body_dm_aa(h1,h2,p2,p1) -= ckl
+             enddo
+             do i = 1, elec_beta_num
+               p2=occ(i,2)
+               h2=p2
+               two_body_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_body_dm_ab(h1,h2,p1,p2) += ckl
+             enddo
+             do i = 1, elec_beta_num
+               p2=occ(i,2)
+               h2=p2
+               two_body_dm_bb(h1,h2,p1,p2) += ckl
+               two_body_dm_bb(h1,h2,p2,p1) -= ckl
+             enddo
+         end select
+       endif
+     enddo
+   enddo
+   ! OMP END DO
+   ! OMP END PARALLEL
+   
+END_PROVIDER
+ 
+ 
diff --git a/plugins/DensityMatrix/det_num.irp.f b/plugins/DensityMatrix/det_num.irp.f
deleted file mode 100644
index 6ba0c1b4..00000000
--- a/plugins/DensityMatrix/det_num.irp.f
+++ /dev/null
@@ -1,56 +0,0 @@
- use bitmasks
-
- BEGIN_PROVIDER [integer, det_num]
-  det_num = 10
- END_PROVIDER
-
- BEGIN_PROVIDER [ integer(bit_kind), det_provider, (N_int,2,det_num)]
-&BEGIN_PROVIDER [ double precision , det_coef_provider, (det_num) ]
- use bitmasks
- implicit none
- integer :: i
- det_provider = 0
- det_provider(1,1,1 ) = #001f ! 0000 0000 0001 1111
- det_provider(1,1,2 ) = #003b ! 0000 0000 0011 1011
- det_provider(1,1,3 ) = #008f ! 0000 0000 1000 1111
- det_provider(1,1,4 ) = #0057 ! 0000 0000 0101 0111
- det_provider(1,1,5 ) = #100f ! 0001 0000 0000 1111
- det_provider(1,1,6 ) = #001f ! 0000 0000 0001 1111
- det_provider(1,1,7 ) = #003b ! 0000 0000 0011 1011
- det_provider(1,1,8 ) = #00c7 ! 0000 0000 1100 0111
- det_provider(1,1,9 ) = #00ab ! 0000 0000 1010 1011
- det_provider(1,1,10) = #0073 ! 0000 0000 0111 0011
- det_provider(1,2,1 ) = #0007 ! 0000 0000 0001 0111
- det_provider(1,2,2 ) = #0023 ! 0000 0000 0010 0011
- det_provider(1,2,3 ) = #0023 ! 0000 0000 0010 0011
- det_provider(1,2,4 ) = #0023 ! 0000 0000 0010 0011
- det_provider(1,2,5 ) = #0015 ! 0000 0000 0001 0101
- det_provider(1,2,6 ) = #000d ! 0000 0000 0000 1101
- det_provider(1,2,7 ) = #0007 ! 0000 0000 0000 0111
- det_provider(1,2,8 ) = #0007 ! 0000 0000 0000 0111
- det_provider(1,2,9 ) = #0007 ! 0000 0000 0000 0111
- det_provider(1,2,10) = #0007 ! 0000 0000 0000 0111
- det_coef_provider = (/   &
-   0.993536117982429D+00, &
-  -0.556089064313864D-01, &
-   0.403074722590178D-01, &
-   0.403074717461626D-01, &
-  -0.340290975461932D-01, &
-  -0.340290958781670D-01, &
-  -0.333949939765448D-01, &
-   0.333418373363987D-01, &
-  -0.316337211787351D-01, &
-  -0.316337207748718D-01  &
-/)
-
-!do i=1,10
-!  call write_bitstring( 6, det_provider(1,1,i), N_int )
-!enddo
-!print *,  ''
-!do i=1,10
-!  call write_bitstring( 6, det_provider(1,2,i), N_int )
-!enddo
-!print *,  ''
-
-
-END_PROVIDER
diff --git a/plugins/DensityMatrix/print_2rdm.irp.f b/plugins/DensityMatrix/print_2rdm.irp.f
new file mode 100644
index 00000000..b4277938
--- /dev/null
+++ b/plugins/DensityMatrix/print_2rdm.irp.f
@@ -0,0 +1,67 @@
+program pouet
+ implicit none
+ integer :: i,j,k,l
+ double precision, external :: get_mo_bielec_integral
+ read_wf = .True.
+ TOUCH read_wf
+
+  double precision :: e(10)
+  e = 0.d0
+
+ print *,  '1RDM ALPHA'
+ do i=1,mo_tot_num
+  do j=1,mo_tot_num
+   print *,  i, j, one_body_dm_mo_alpha(i,j,1) 
+   e(4) += one_body_dm_mo_alpha(i,j,1) * mo_mono_elec_integral(i,j)
+  enddo
+ enddo
+
+ print *,  '1RDM BETA'
+ do i=1,mo_tot_num
+  do j=1,mo_tot_num
+   print *,  i, j, one_body_dm_mo_beta(i,j,1) 
+   e(4) += one_body_dm_mo_beta(i,j,1) * mo_mono_elec_integral(i,j)
+  enddo
+ enddo
+
+ print *,  '2RDM ALPHA ALPHA'
+ do i=1,mo_tot_num
+   do j=1,mo_tot_num
+     do k=1,mo_tot_num
+       do l=1,mo_tot_num
+         print *,  i, j, k, l, two_body_dm_aa(i,j,k,l)
+         e(1) += two_body_dm_aa(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
+       enddo
+     enddo
+   enddo
+ enddo
+
+ print *,  '2RDM BETA BETA'
+ do i=1,mo_tot_num
+   do j=1,mo_tot_num
+     do k=1,mo_tot_num
+       do l=1,mo_tot_num
+         print *,  i, j, k, l, two_body_dm_bb(i,j,k,l)
+         e(2) += two_body_dm_bb(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
+       enddo
+     enddo
+   enddo
+ enddo
+
+ print *,  '2RDM ALPHA BETA'
+ do i=1,mo_tot_num
+   do j=1,mo_tot_num
+     do k=1,mo_tot_num
+       do l=1,mo_tot_num
+         print *,  i, j, k, l, two_body_dm_ab(i,j,k,l)
+         e(3) += two_body_dm_ab(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
+       enddo
+     enddo
+   enddo
+ enddo
+
+ print *,  ''
+ print *,  'Energy ', sum(e(1:4)) + nuclear_repulsion
+
+
+end