diff --git a/devel/cassdtq/EZFIO.cfg b/devel/cassdtq/EZFIO.cfg
new file mode 100644
index 0000000..e2af4e5
--- /dev/null
+++ b/devel/cassdtq/EZFIO.cfg
@@ -0,0 +1,25 @@
+[energy]
+type: double precision
+doc: Calculated Selected CASSD energy
+interface: ezfio
+size: (determinants.n_states)
+
+[energy_pt2]
+type: double precision
+doc: Calculated CASSD energy + PT2
+interface: ezfio
+size: (determinants.n_states)
+
+
+[do_ddci]
+type: logical
+doc: If true, remove purely inactive double excitations
+interface: ezfio,provider,ocaml
+default: False
+
+[do_only_1h1p]
+type: logical
+doc: If true, do only one hole/one particle excitations
+interface: ezfio,provider,ocaml
+default: False
+
diff --git a/devel/cassdtq/NEED b/devel/cassdtq/NEED
new file mode 100644
index 0000000..9afaa86
--- /dev/null
+++ b/devel/cassdtq/NEED
@@ -0,0 +1,3 @@
+cipsi                                                                                         
+selectors_full
+
diff --git a/devel/cassdtq/README.rst b/devel/cassdtq/README.rst
new file mode 100644
index 0000000..8d3025d
--- /dev/null
+++ b/devel/cassdtq/README.rst
@@ -0,0 +1,5 @@
+======
+cisdtq
+======
+
+Selected CI in the CISDTQ space.
diff --git a/devel/cassdtq/cassdtq.irp.f b/devel/cassdtq/cassdtq.irp.f
new file mode 100644
index 0000000..680720e
--- /dev/null
+++ b/devel/cassdtq/cassdtq.irp.f
@@ -0,0 +1,22 @@
+program cassdtq
+  implicit none
+  BEGIN_DOC
+  ! Selected CISDTQ with stochastic selection and PT2.
+  END_DOC
+
+  if (.not.is_zmq_slave) then
+    PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
+
+    if (do_pt2) then
+      call run_stochastic_cipsi
+    else
+      call run_cipsi
+    endif
+
+  else
+    PROVIDE mo_two_e_integrals_in_map
+
+    call run_slave_cipsi
+
+  endif
+end
diff --git a/devel/cassdtq/class.irp.f b/devel/cassdtq/class.irp.f
new file mode 100644
index 0000000..6d1838a
--- /dev/null
+++ b/devel/cassdtq/class.irp.f
@@ -0,0 +1,8 @@
+BEGIN_PROVIDER [ logical, do_only_cas  ]
+ implicit none
+ BEGIN_DOC
+ ! In the CASSDTQ case, all those are always false
+ END_DOC
+ do_only_cas  = .False.
+END_PROVIDER
+
diff --git a/devel/cassdtq/generators.irp.f b/devel/cassdtq/generators.irp.f
new file mode 100644
index 0000000..66ffe10
--- /dev/null
+++ b/devel/cassdtq/generators.irp.f
@@ -0,0 +1,78 @@
+use bitmasks
+
+logical function is_a_generator(det)
+  use bitmasks
+  implicit none
+  integer(bit_kind), intent(in)  :: det(N_int,2)
+  integer :: na, nb
+  integer, external :: number_of_holes, number_of_particles
+  provide N_int
+  na = number_of_holes(det)
+  nb = number_of_particles(det)
+  is_a_generator = (na <= 2) .and. (nb <= 2)
+end
+
+BEGIN_PROVIDER [ integer, N_det_generators ]
+  implicit none
+  BEGIN_DOC
+  ! Number of generator detetrminants
+  END_DOC
+  integer                        :: i,k,l
+  logical, external              :: is_a_generator
+  provide N_int
+  call write_time(6)
+  N_det_generators = 0
+  do i=1,N_det
+    if (is_a_generator(psi_det_sorted(1,1,i))) then
+      N_det_generators += 1
+    endif
+  enddo
+  N_det_generators = max(N_det_generators,1)
+  call write_int(6,N_det_generators,'Number of generators')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order, (psi_det_size) ]
+  implicit none
+  BEGIN_DOC
+  ! For Single reference wave functions, the generator is the
+  ! Hartree-Fock determinant
+  END_DOC
+  integer                        :: i, k, l, m
+  logical, external              :: is_a_generator
+  integer, allocatable           :: nongen(:)
+  integer                        :: inongen
+
+  allocate(nongen(N_det))
+
+  inongen = 0
+  m=0
+  do i=1,N_det
+    if (is_a_generator(psi_det_sorted(1,1,i))) then
+      m = m+1
+      psi_det_sorted_gen_order(i) = m
+      do k=1,N_int
+        psi_det_generators(k,1,m) = psi_det_sorted(k,1,i)
+        psi_det_generators(k,2,m) = psi_det_sorted(k,2,i)
+      enddo
+      psi_coef_generators(m,:) = psi_coef_sorted(i,:)
+    else
+      inongen += 1
+      nongen(inongen) = i
+    endif
+  enddo
+
+  psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators)
+  psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :)
+  do i=1,inongen
+    psi_det_sorted_gen_order(nongen(i)) = N_det_generators+i
+    psi_det_sorted_gen(:,:,N_det_generators+i) = psi_det_sorted(:,:,nongen(i))
+    psi_coef_sorted_gen(N_det_generators+i, :) = psi_coef_sorted(nongen(i),:)
+  end do
+
+END_PROVIDER
+
+
diff --git a/devel/cassdtq/save_energy.irp.f b/devel/cassdtq/save_energy.irp.f
new file mode 100644
index 0000000..b0441eb
--- /dev/null
+++ b/devel/cassdtq/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_cassdtq_energy(E(1:N_states))
+  call ezfio_set_cassdtq_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end