diff --git a/devel/mpn/.gitignore b/devel/mpn/.gitignore
new file mode 100644
index 0000000..1561915
--- /dev/null
+++ b/devel/mpn/.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/mpn/EZFIO.cfg b/devel/mpn/EZFIO.cfg
new file mode 100644
index 0000000..50e5532
--- /dev/null
+++ b/devel/mpn/EZFIO.cfg
@@ -0,0 +1,12 @@
+[energy]
+type: double precision
+doc: Calculated Selected |FCI| energy
+interface: ezfio
+size: (determinants.n_states)
+
+[mp_order]
+type: integer
+doc: Max order of MPn
+interface: ezfio, provider, ocaml
+default: 4
+
diff --git a/devel/mpn/NEED b/devel/mpn/NEED
new file mode 100644
index 0000000..b2918cc
--- /dev/null
+++ b/devel/mpn/NEED
@@ -0,0 +1,3 @@
+davidson_undressed
+hartree_fock
+determinants
diff --git a/devel/mpn/README.rst b/devel/mpn/README.rst
new file mode 100644
index 0000000..4967020
--- /dev/null
+++ b/devel/mpn/README.rst
@@ -0,0 +1,4 @@
+===
+mpn
+===
+
diff --git a/devel/mpn/energies.irp.f b/devel/mpn/energies.irp.f
new file mode 100644
index 0000000..c3a1448
--- /dev/null
+++ b/devel/mpn/energies.irp.f
@@ -0,0 +1,23 @@
+BEGIN_PROVIDER [ double precision, energy_det_i, (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Fock Energy of determinant |I> (sum of epsilon_i) 
+ END_DOC
+ integer :: i, k, n
+ integer :: list(elec_alpha_num)
+
+ do k=1,N_det
+   call bitstring_to_list(psi_det(1,1,k), list, n, N_int)
+   energy_det_i(k) = 0.d0
+   do i=1,n
+     energy_det_i(k) += fock_matrix_diag_mo(list(i))
+   enddo
+   call bitstring_to_list(psi_det(1,2,k), list, n, N_int)
+   do i=1,n
+     energy_det_i(k) += fock_matrix_diag_mo(list(i))
+   enddo
+ enddo
+END_PROVIDER
+
+
+
diff --git a/devel/mpn/mpn.irp.f b/devel/mpn/mpn.irp.f
new file mode 100644
index 0000000..a988913
--- /dev/null
+++ b/devel/mpn/mpn.irp.f
@@ -0,0 +1,101 @@
+program mpn
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  integer :: i, k, l
+  double precision, allocatable :: c_pert(:,:)
+  double precision, allocatable :: e_pert(:)
+  double precision, allocatable :: hc(:), s2(:)
+
+  n_states_diag = 1
+  TOUCH n_states_diag
+  call generate_fci_space
+  allocate(c_pert(N_det,0:mp_order))
+  allocate(s2(N_det))
+  allocate(e_pert(mp_order+1))
+  e_pert = 0.d0
+  c_pert(:,:) = 0.d0
+  c_pert(1,0) = 1.d0
+
+double precision :: hij
+
+  do k=1,mp_order
+    ! H_ij C^(k-1)
+    call h_s2_u_0_nstates_zmq(c_pert(1,k),s2,c_pert(1,k-1),1,N_det)
+    e_pert(k) = c_pert(1,k)
+    print *, k, e_pert(k), sum(e_pert) + nuclear_repulsion
+
+    c_pert(1,k) = 0.d0
+    c_pert(:,k) = -c_pert(:,k)
+    do l=1,k-1
+      do i=2,N_det
+        c_pert(i,k) = c_pert(i,k) + e_pert(l) * c_pert(i,k-l)
+      enddo
+    enddo
+    do i=2,N_det
+      c_pert(i,k) = c_pert(i,k) + energy_det_i(i) * c_pert(i,k-1)
+    enddo
+    do i=2,N_det
+      c_pert(i,k) = c_pert(i,k) / (energy_det_i(i) - energy_det_i(1))
+    enddo
+
+  enddo
+
+
+end
+
+subroutine generate_fci_space
+  use bitmasks
+  implicit none
+  integer :: i, sze
+  integer(bit_kind) :: o(N_int,2)
+
+  if (mo_num > 64) then
+     stop 'No more than 64 MOs'
+  endif
+  o(:,1) = full_ijkl_bitmask(:)
+  o(:,2) = 0_bit_kind
+
+  call configuration_to_dets_size(o,n_det_alpha_unique,elec_alpha_num,N_int)
+  TOUCH n_det_alpha_unique
+
+  integer :: k,n,m, t, t1, t2
+  integer(bit_kind) :: u
+  k=0
+  n = elec_alpha_num
+  m = mo_num - n
+  u = shiftl(1_bit_kind,n) -1
+  do while (u < shiftl(1_bit_kind,n+m))
+    k = k+1
+    psi_det_alpha_unique(1,k) = u
+    t = ior(u,u-1)
+    t1 = t+1
+    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
+    u = ior(t1,t2)
+  enddo
+
+
+  call configuration_to_dets_size(o,n_det_beta_unique,elec_beta_num,N_int)
+  TOUCH n_det_beta_unique
+
+  k=0
+  n = elec_beta_num
+  m = mo_num - n
+  u = shiftl(1_bit_kind,n) -1
+  do while (u < shiftl(1_bit_kind,n+m))
+    k = k+1
+    psi_det_beta_unique(1,k) = u
+    t = ior(u,u-1)
+    t1 = t+1
+    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
+    u = ior(t1,t2)
+  enddo
+
+  call generate_all_alpha_beta_det_products
+
+  print *, N_det
+
+
+end
+
diff --git a/devel/svdwf/svdwf.irp.f b/devel/svdwf/svdwf.irp.f
index dfdcc80..1f67022 100644
--- a/devel/svdwf/svdwf.irp.f
+++ b/devel/svdwf/svdwf.irp.f
@@ -19,16 +19,27 @@ subroutine run
             Vt(n_det_beta_unique, n_det_beta_unique),  &
             D(max(n_det_beta_unique,n_det_alpha_unique)) )
 
-  A = 0.D0
+  do j=1,n_det_beta_unique
+    do i=1,n_det_alpha_unique
+      A(i,j)  = 0.D0
+    enddo
+  enddo
+
   do k=1,N_det
     i = psi_bilinear_matrix_rows(k)
     j = psi_bilinear_matrix_columns(k)
     A(i,j) = psi_bilinear_matrix_values(k,1)
   enddo
 
-  call randomized_svd(A, size(A,1), &
+  if (N_det == 1) then
+    D(1) = 1.d0
+    U(1,1) = 1.d0
+    Vt(1,1) = 1.d0
+  else
+    call randomized_svd(A, size(A,1), &
            U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
-           6,1000)
+           6,min(n_det_beta_unique,1000))
+  endif
 
   entropy = 0.d0
   do i=1,n_det_beta_unique
@@ -41,6 +52,7 @@ subroutine run
   enddo
   print *, 'threshold: ', 2.858 * D(k/2)
   print *, 'Entropy  : ', entropy
+
 !  do i=1,n_det_alpha_unique
 !    print '(I6,4(X,F12.8))', i, U(i,1:4)
 !  enddo
@@ -48,4 +60,5 @@ subroutine run
 !  do i=1,n_det_beta_unique
 !    print '(I6,4(X,F12.8))', i, Vt(1:4,i)
 !  enddo
+
 end