diff --git a/stable/mp2/H_apply.irp.f b/stable/mp2/H_apply.irp.f
new file mode 100644
index 0000000..4d6e5a4
--- /dev/null
+++ b/stable/mp2/H_apply.irp.f
@@ -0,0 +1,15 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python2 ]
+from generate_h_apply import *
+from perturbation import perturbations
+
+s = H_apply("mp2")
+s.set_perturbation("Moller_plesset")
+#s.set_perturbation("epstein_nesbet")
+print s
+
+s = H_apply("mp2_selection")
+s.set_selection_pt2("Moller_Plesset")
+print s
+END_SHELL
+
diff --git a/stable/mp2/NEED b/stable/mp2/NEED
new file mode 100644
index 0000000..6eaf5b9
--- /dev/null
+++ b/stable/mp2/NEED
@@ -0,0 +1,6 @@
+generators_full
+selectors_full
+determinants
+davidson
+davidson_undressed
+perturbation
diff --git a/stable/mp2/README.rst b/stable/mp2/README.rst
new file mode 100644
index 0000000..192a75f
--- /dev/null
+++ b/stable/mp2/README.rst
@@ -0,0 +1,4 @@
+===
+mp2
+===
+
diff --git a/stable/mp2/mp2.irp.f b/stable/mp2/mp2.irp.f
new file mode 100644
index 0000000..b8e0cc4
--- /dev/null
+++ b/stable/mp2/mp2.irp.f
@@ -0,0 +1,21 @@
+program mp2
+  call run
+end
+
+subroutine run
+  implicit none
+  double precision, allocatable  :: pt2(:), norm_pert(:)
+  double precision               :: H_pert_diag, E_old
+  integer                        :: N_st, iter
+  PROVIDE Fock_matrix_diag_mo H_apply_buffer_allocated
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st))
+  E_old = HF_energy
+  call H_apply_mp2(pt2, norm_pert, H_pert_diag,  N_st)
+  print *,  'N_det    = ', N_det
+  print *,  'N_states = ', N_states
+  print *,  'MP2      = ', pt2
+  print *,  'E        = ', E_old
+  print *,  'E+MP2    = ', E_old+pt2
+  deallocate(pt2,norm_pert)
+end