diff --git a/devel/density_matrix/.gitignore b/devel/density_matrix/.gitignore
new file mode 100644
index 0000000..7ac9fbf
--- /dev/null
+++ b/devel/density_matrix/.gitignore
@@ -0,0 +1,5 @@
+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
new file mode 100644
index 0000000..d3d4d2c
--- /dev/null
+++ b/devel/density_matrix/NEED
@@ -0,0 +1 @@
+determinants
diff --git a/devel/density_matrix/README.rst b/devel/density_matrix/README.rst
new file mode 100644
index 0000000..ded0613
--- /dev/null
+++ b/devel/density_matrix/README.rst
@@ -0,0 +1,5 @@
+====================
+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
new file mode 100644
index 0000000..51ab45c
--- /dev/null
+++ b/devel/density_matrix/density_matrix_array.irp.f
@@ -0,0 +1,115 @@
+ 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
new file mode 100644
index 0000000..682222a
--- /dev/null
+++ b/devel/density_matrix/print_2rdm.irp.f
@@ -0,0 +1,43 @@
+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
new file mode 100644
index 0000000..1b94fe5
--- /dev/null
+++ b/devel/density_matrix/print_2rdm_decomposed.irp.f
@@ -0,0 +1,78 @@
+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/stable/champ/.gitignore b/stable/champ/.gitignore
new file mode 100644
index 0000000..7ac9fbf
--- /dev/null
+++ b/stable/champ/.gitignore
@@ -0,0 +1,5 @@
+IRPF90_temp/
+IRPF90_man/
+irpf90.make
+irpf90_entities
+tags
\ No newline at end of file
diff --git a/stable/mp2/.gitignore b/stable/mp2/.gitignore
new file mode 100644
index 0000000..1561915
--- /dev/null
+++ b/stable/mp2/.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