diff --git a/src/cisd/EZFIO.cfg b/src/cisd/EZFIO.cfg
index 4565d2df..688f802a 100644
--- a/src/cisd/EZFIO.cfg
+++ b/src/cisd/EZFIO.cfg
@@ -5,3 +5,11 @@ interface: ezfio
 size: (determinants.n_states)
 
 
+
+[lcc_energy]
+type: double precision
+doc: lccsd energy 
+interface: ezfio
+size: (determinants.n_states)
+
+
diff --git a/src/cisd/NEED b/src/cisd/NEED
index d9ad3efc..616d021e 100644
--- a/src/cisd/NEED
+++ b/src/cisd/NEED
@@ -1,3 +1,4 @@
 selectors_full
 single_ref_method
 davidson_undressed
+dav_general_mat
diff --git a/src/cisd/lccsd.irp.f b/src/cisd/lccsd.irp.f
new file mode 100644
index 00000000..919c5aaa
--- /dev/null
+++ b/src/cisd/lccsd.irp.f
@@ -0,0 +1,95 @@
+program lccsd
+  implicit none
+  BEGIN_DOC
+! Linerarized CCSD 
+!
+  ! This program takes a reference Slater determinant of ROHF-like occupancy,
+  !
+  ! 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`).
+  !
+  ! This program can be useful in many cases:
+  !
+  ! * **Ground state calculation**: if even after a :c:func:`cis` calculation, natural
+  !   orbitals (see :c:func:`save_natorb`) and then :c:func:`scf` optimization, you are not sure to have the lowest scf
+  !   solution,
+  !   do the same strategy with the :c:func:`cisd` executable instead of the :c:func:`cis` exectuable to generate the natural
+  !   orbitals as a guess for the :c:func:`scf`.
+  !
+  !
+  !
+  ! * **Excited states calculations**: the lowest excited states are much likely to
+  !   be dominanted by single- or double-excitations.
+  !   Therefore, running a :c:func:`cisd` will save the "n_states" lowest states within
+  !   the CISD space
+  !   in the |EZFIO| directory, which can afterward be used as guess wave functions
+  !   for a further multi-state fci calculation if you specify "read_wf" = True
+  !   before running the fci executable (see :option:`determinants read_wf`).
+  !   Also, if you specify "s2_eig" = True, the cisd will only retain states
+  !   having the good value :math:`S^2` value
+  !   (see :option:`determinants expected_s2` and :option:`determinants s2_eig`).
+  !   If "s2_eig" = False, it will take the lowest n_states, whatever
+  !   multiplicity they are.
+  !
+  !
+  !
+  !   Note: if you would like to discard some orbitals, use
+  !   :ref:`qp_set_mo_class` to specify:
+  !
+  !   * "core" orbitals which will be always doubly occupied
+  !
+  !   * "act" orbitals where an electron can be either excited from or to
+  !
+  !   * "del" orbitals which will be never occupied
+  !
+  END_DOC
+  PROVIDE N_states
+  read_wf = .False.
+  TOUCH read_wf
+  call run
+end
+
+subroutine run
+  implicit none
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  call get_lccsd_2
+end
+
+subroutine get_lccsd_2
+  implicit none
+  integer                        :: i,k
+  double precision               :: cisdq(N_states), delta_e
+  double precision,external      :: diag_h_mat_elem
+  psi_coef = lccsd_coef
+  SOFT_TOUCH psi_coef
+  call save_wavefunction_truncated(save_threshold)
+  call ezfio_set_cisd_lcc_energy(lccsd_energies)
+
+  print *,  'N_det = ', N_det
+  print*,''
+  print*,'******************************'
+  print *,  'LCCSD Energies'
+  do i = 1,N_states
+    print *,  i, lccsd_energies(i)
+  enddo
+  if (N_states > 1) then
+      print*,'******************************'
+      print*,'Excitation energies (au)    (LCCSD)'
+      do i = 2, N_states
+        print*, i ,lccsd_energies(i) - lccsd_energies(1)
+      enddo
+      print*,''
+      print*,'******************************'
+      print*,'Excitation energies (eV)    (LCCSD)'
+      do i = 2, N_states
+        print*, i ,(lccsd_energies(i) - lccsd_energies(1)) * ha_to_ev
+      enddo
+  endif
+
+end
diff --git a/src/cisd/lccsd_prov.irp.f b/src/cisd/lccsd_prov.irp.f
new file mode 100644
index 00000000..b071a8f8
--- /dev/null
+++ b/src/cisd/lccsd_prov.irp.f
@@ -0,0 +1,44 @@
+ BEGIN_PROVIDER [ double precision, lccsd_coef, (N_det, N_states)]
+&BEGIN_PROVIDER [ double precision, lccsd_energies, (N_states)]
+ implicit none
+ double precision, allocatable :: Dress_jj(:), H_jj(:), u_in(:,:)
+ double precision :: ebefore, eafter, ecorr, thresh
+ integer :: i,it
+ logical :: converged
+ external H_u_0_nstates_openmp
+ allocate(Dress_jj(N_det),H_jj(N_det),u_in(N_det,N_states_diag))
+ thresh = 1.d-6
+ converged = .False.
+ Dress_jj = 0.d0
+ u_in = 0.d0
+ it = 0
+ ! initial guess
+ do i = 1, N_states_diag
+  u_in(i,i) = 1.d0
+ enddo
+ do i = 1,N_det
+  call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,H_jj(i))
+ enddo
+ ebefore = H_jj(1)
+ do while (.not.converged)
+  it += 1
+  print*,'N_det = ',N_det
+  call davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,lccsd_energies,& 
+                   N_det,N_states,N_states_diag,converged,H_u_0_nstates_openmp)
+  ecorr = lccsd_energies(1) - H_jj(1) 
+  print*,'---------------------'
+  print*,'it = ',it
+  print*,'ecorr = ',ecorr
+  Dress_jj(1) = 0.d0
+  do i = 2, N_det
+    Dress_jj(i) = ecorr
+  enddo
+  eafter = lccsd_energies(1)
+  converged = (dabs(eafter - ebefore).lt.thresh)
+  ebefore = eafter
+ enddo
+ do i = 1, N_states
+  lccsd_coef(1:N_det,i) = u_in(1:N_det,i)
+ enddo
+
+END_PROVIDER