diff --git a/devel/general_mrci/.gitignore b/devel/general_mrci/.gitignore
new file mode 100644
index 0000000..ed67ed1
--- /dev/null
+++ b/devel/general_mrci/.gitignore
@@ -0,0 +1,6 @@
+IRPF90_temp/
+IRPF90_man/
+irpf90.make
+irpf90_entities
+tags
+mrci
diff --git a/devel/general_mrci/EZFIO.cfg b/devel/general_mrci/EZFIO.cfg
new file mode 100644
index 0000000..ddd57ef
--- /dev/null
+++ b/devel/general_mrci/EZFIO.cfg
@@ -0,0 +1,23 @@
+[energy]
+type: double precision
+doc: Calculated Selected |FCI| energy
+interface: ezfio
+size: (determinants.n_states)
+
+[energy_pt2]
+type: double precision
+doc: Calculated |FCI| energy + |PT2|
+interface: ezfio
+size: (determinants.n_states)
+
+[do_singles]
+type: logical
+doc: If true, add the singly excited determinants
+interface: ezfio, provider, ocaml
+default: true
+
+[do_doubles]
+type: logical
+doc: If true, add the doubly excited determinants
+interface: ezfio, provider, ocaml
+default: true
diff --git a/devel/general_mrci/NEED b/devel/general_mrci/NEED
new file mode 100644
index 0000000..42992ac
--- /dev/null
+++ b/devel/general_mrci/NEED
@@ -0,0 +1,3 @@
+selectors_full
+generators_full
+davidson_undressed
diff --git a/devel/general_mrci/README.rst b/devel/general_mrci/README.rst
new file mode 100644
index 0000000..d3d9a62
--- /dev/null
+++ b/devel/general_mrci/README.rst
@@ -0,0 +1,6 @@
+============
+general_mrci
+============
+
+Allows the computation of singles and doubles on top of an arbitrary reference wave function.
+
diff --git a/devel/general_mrci/class.irp.f b/devel/general_mrci/class.irp.f
new file mode 100644
index 0000000..b4a68ac
--- /dev/null
+++ b/devel/general_mrci/class.irp.f
@@ -0,0 +1,12 @@
+ BEGIN_PROVIDER [ logical, do_only_1h1p ]
+&BEGIN_PROVIDER [ logical, do_only_cas  ]
+&BEGIN_PROVIDER [ logical, do_ddci ]
+ implicit none
+ BEGIN_DOC
+ ! In the FCI case, all those are always false
+ END_DOC
+ do_only_1h1p = .False.
+ do_only_cas  = .False.
+ do_ddci = .False.
+END_PROVIDER
+
diff --git a/devel/general_mrci/h_apply.irp.f b/devel/general_mrci/h_apply.irp.f
new file mode 100644
index 0000000..bb0d944
--- /dev/null
+++ b/devel/general_mrci/h_apply.irp.f
@@ -0,0 +1,16 @@
+! Generates subroutine H_apply_cisd
+! ----------------------------------
+
+BEGIN_SHELL [ /usr/bin/env python3 ]
+from generate_h_apply import H_apply
+H = H_apply("cisd")
+print(H)
+
+H = H_apply("cis",do_double_exc=False)
+print(H)
+
+H = H_apply("cid",do_mono_exc=False)
+print(H)
+
+END_SHELL
+
diff --git a/devel/general_mrci/mrci.irp.f b/devel/general_mrci/mrci.irp.f
new file mode 100644
index 0000000..8b56f25
--- /dev/null
+++ b/devel/general_mrci/mrci.irp.f
@@ -0,0 +1,77 @@
+program cisd
+  implicit none
+  BEGIN_DOC
+  ! Multi-Reference Configuration Interaction with Single and Double excitations.
+  !
+  ! This program takes a reference wave function as input and
+  ! and performs all single and double excitations on top of it, disregarding
+  ! spatial symmetry and compute the "n_states" lowest eigenstates of that CI
+  ! matrix (see :option:`determinants n_states`).
+  !
+  END_DOC
+  PROVIDE N_states
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call run
+end
+
+subroutine run
+  implicit none
+  integer                        :: i,k
+  double precision               :: cisdq(N_states), delta_e
+  double precision,external      :: diag_h_mat_elem
+
+  if(do_singles)then
+    if(do_doubles)then
+      call H_apply_cisd
+    else
+      call H_apply_cis
+    endif
+  else
+    if(do_doubles)then
+      call H_apply_cid
+    else
+      stop "do_singles and do_doubles are both False"
+    endif
+  endif
+
+  psi_coef = ci_eigenvectors
+  SOFT_TOUCH psi_coef
+  call save_wavefunction
+  call ezfio_set_general_mrci_energy(CI_energy)
+
+  do i = 1,N_states
+    k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
+    delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
+    cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
+  enddo
+  print *,  'N_det = ', N_det
+  print*,''
+  print*,'******************************'
+  print *,  'MR-CI Energies'
+  do i = 1,N_states
+    print *,  i, CI_energy(i)
+  enddo
+  print*,''
+  print*,'******************************'
+  print *,  'MR-CI+Q Energies'
+  do i = 1,N_states
+    print *,  i, cisdq(i)
+  enddo
+  if (N_states > 1) then
+    print*,''
+    print*,'******************************'
+    print*,'Excitation energies (au)    (CISD+Q)'
+    do i = 2, N_states
+      print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1)
+    enddo
+    print*,''
+    print*,'******************************'
+    print*,'Excitation energies (eV)    (CISD+Q)'
+    do i = 2, N_states
+      print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, &
+        (cisdq(i) - cisdq(1)) / 0.0367502d0
+    enddo
+  endif
+
+end
diff --git a/devel/general_mrci/save_energy.irp.f b/devel/general_mrci/save_energy.irp.f
new file mode 100644
index 0000000..b553f41
--- /dev/null
+++ b/devel/general_mrci/save_energy.irp.f
@@ -0,0 +1,9 @@
+subroutine save_energy(E,pt2)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E(N_states), pt2(N_states)
+  call ezfio_set_general_mrci_energy(E(1:N_states))
+  call ezfio_set_general_mrci_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end