Program FOBO-SCF works

plugins/FOBOCI/EZFIO.cfg

@@ -1,6 +1,13 @@
-[threshold_singles]
+[threshold_lmct]
 type: double precision
-doc: threshold to select the pertinent single excitations at second order
+doc: threshold to select the pertinent LMCT excitations at second order
+interface: ezfio,provider,ocaml
+default: 0.01
+
+
+[threshold_mlct]
+type: double precision
+doc: threshold to select the pertinent MLCT excitations at second order
 interface: ezfio,provider,ocaml
 default: 0.01
 
@@ -16,6 +23,20 @@ doc: if true, you do the FOBOCI calculation perturbatively
 interface: ezfio,provider,ocaml
 default: .False.
 
+
+[speed_up_convergence_foboscf]
+type: logical
+doc: if true, the threshold of the FOBO-SCF algorithms are increased with the iterations
+interface: ezfio,provider,ocaml
+default: .True.
+
+
+[dressing_2h2p]
+type: logical
+doc: if true, you do dress with 2h2p excitations each FOBOCI matrix
+interface: ezfio,provider,ocaml
+default: .False.
+
 [second_order_h]
 type: logical
 doc: if true, you do the FOBOCI calculation using second order intermediate Hamiltonian

plugins/FOBOCI/all_singles.irp.f

@@ -66,7 +66,7 @@ subroutine all_single
    print*,'E        = ',CI_energy(i)
    print*,'S^2      = ',CI_eigenvectors_s2(i)
   enddo
-  do i = 1, 2
+  do i = 1, max(2,N_det_generators)
    print*,'psi_coef = ',psi_coef(i,1)
   enddo
   deallocate(pt2,norm_pert,E_before)

plugins/FOBOCI/collect_all_lmct.irp.f

@@ -92,7 +92,7 @@ subroutine collect_lmct_mlct(hole_particle,n_couples)
   iorb = list_act(i)
   do j = 1, n_inact_orb
    jorb = list_inact(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then
     n_couples +=1
     hole_particle(n_couples,1) = jorb
     hole_particle(n_couples,2) = iorb
@@ -102,7 +102,7 @@ subroutine collect_lmct_mlct(hole_particle,n_couples)
   enddo
   do j = 1, n_virt_orb
    jorb = list_virt(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then
     n_couples +=1
     hole_particle(n_couples,1) = iorb
     hole_particle(n_couples,2) = jorb

plugins/FOBOCI/corr_energy_2h2p.irp.f

@@ -1,8 +1,16 @@
  BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h2p_ab_2_orb, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h2p_bb_2_orb, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_a, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_b, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_double, (mo_tot_num,mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)]
 &BEGIN_PROVIDER [ double precision, total_corr_e_2h2p]
  use bitmasks
+ print*,''
+ print*,'Providing the 2h2p correlation energy'
+ print*,''
  implicit none
  integer(bit_kind) :: key_tmp(N_int,2)
  integer :: i,j,k,l
@@ -12,9 +20,14 @@
  integer :: i_ok,ispin
  ! Alpha - Beta correlation energy
  total_corr_e_2h2p = 0.d0
+ corr_energy_2h2p_ab_2_orb = 0.d0
+ corr_energy_2h2p_bb_2_orb = 0.d0
  corr_energy_2h2p_per_orb_ab = 0.d0
  corr_energy_2h2p_per_orb_aa = 0.d0
  corr_energy_2h2p_per_orb_bb = 0.d0
+ corr_energy_2h2p_for_1h1p_a = 0.d0
+ corr_energy_2h2p_for_1h1p_b = 0.d0
+ corr_energy_2h2p_for_1h1p_double = 0.d0
  do i = 1, n_inact_orb    ! beta 
   i_hole = list_inact(i)
   do k = 1, n_virt_orb    ! beta
@@ -36,8 +49,24 @@
      hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
      contrib = hij*hij/delta_e
      total_corr_e_2h2p += contrib
+     ! Single orbital contribution
      corr_energy_2h2p_per_orb_ab(i_hole) += contrib
      corr_energy_2h2p_per_orb_ab(k_part) += contrib
+     ! Couple of orbital contribution for the single 1h1p
+     corr_energy_2h2p_for_1h1p_a(j_hole,l_part) += contrib
+     corr_energy_2h2p_for_1h1p_a(l_part,j_hole) += contrib
+     corr_energy_2h2p_for_1h1p_b(j_hole,l_part) += contrib
+     corr_energy_2h2p_for_1h1p_b(l_part,j_hole) += contrib
+     ! Couple of orbital contribution for the double 1h1p
+     corr_energy_2h2p_for_1h1p_double(i_hole,l_part) += contrib
+     corr_energy_2h2p_for_1h1p_double(l_part,i_hole) += contrib
+
+     corr_energy_2h2p_ab_2_orb(i_hole,j_hole) += contrib
+     corr_energy_2h2p_ab_2_orb(j_hole,i_hole) += contrib
+     corr_energy_2h2p_ab_2_orb(i_hole,k_part) += contrib
+     corr_energy_2h2p_ab_2_orb(k_part,i_hole) += contrib
+     corr_energy_2h2p_ab_2_orb(k_part,l_part) += contrib
+     corr_energy_2h2p_ab_2_orb(l_part,k_part) += contrib
     enddo 
    enddo
   enddo
@@ -65,8 +94,12 @@
      hij = hij  - exc 
      contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
      total_corr_e_2h2p += contrib
+     ! Single orbital contribution
      corr_energy_2h2p_per_orb_aa(i_hole) += contrib
      corr_energy_2h2p_per_orb_aa(k_part) += contrib
+     ! Couple of orbital contribution for the single 1h1p
+     corr_energy_2h2p_for_1h1p_a(i_hole,k_part) += contrib
+     corr_energy_2h2p_for_1h1p_a(k_part,i_hole) += contrib
     enddo 
    enddo
   enddo
@@ -94,8 +127,20 @@
      hij = hij  - exc 
      contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
      total_corr_e_2h2p += contrib
+     ! Single orbital contribution
      corr_energy_2h2p_per_orb_bb(i_hole) += contrib
      corr_energy_2h2p_per_orb_bb(k_part) += contrib
+     corr_energy_2h2p_for_1h1p_b(i_hole,k_part) += contrib
+     corr_energy_2h2p_for_1h1p_b(k_part,i_hole) += contrib
+
+     ! Two particle correlation energy
+     corr_energy_2h2p_bb_2_orb(i_hole,j_hole) += contrib
+     corr_energy_2h2p_bb_2_orb(j_hole,i_hole) += contrib
+     corr_energy_2h2p_bb_2_orb(i_hole,k_part) += contrib
+     corr_energy_2h2p_bb_2_orb(k_part,i_hole) += contrib
+     corr_energy_2h2p_bb_2_orb(k_part,l_part) += contrib
+     corr_energy_2h2p_bb_2_orb(l_part,k_part) += contrib
+
     enddo 
    enddo
   enddo
@@ -103,7 +148,11 @@
 
 END_PROVIDER 
 
- BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)]
+ BEGIN_PROVIDER [double precision, corr_energy_2h1p_ab_bb_per_2_orb, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_a, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_b, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_double, (mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)]
 &BEGIN_PROVIDER [ double precision, total_corr_e_2h1p]
@@ -120,6 +169,10 @@ END_PROVIDER
  corr_energy_2h1p_per_orb_ab = 0.d0
  corr_energy_2h1p_per_orb_aa = 0.d0
  corr_energy_2h1p_per_orb_bb = 0.d0
+ corr_energy_2h1p_ab_bb_per_2_orb = 0.d0
+ corr_energy_2h1p_for_1h1p_a = 0.d0
+ corr_energy_2h1p_for_1h1p_b = 0.d0
+ corr_energy_2h1p_for_1h1p_double = 0.d0
  do i = 1, n_inact_orb
   i_hole = list_inact(i)
   do k = 1, n_act_orb
@@ -141,6 +194,7 @@ END_PROVIDER
      hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
      contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
      total_corr_e_2h1p += contrib
+     corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib
      corr_energy_2h1p_per_orb_ab(i_hole) += contrib
      corr_energy_2h1p_per_orb_ab(l_part) += contrib
     enddo 
@@ -199,6 +253,7 @@ END_PROVIDER
      delta_e  = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
      hij = hij  - exc 
      contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
+     corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib
 
      total_corr_e_2h1p += contrib
      corr_energy_2h1p_per_orb_bb(i_hole) += contrib
@@ -212,6 +267,7 @@ END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)]
+&BEGIN_PROVIDER [double precision, corr_energy_1h2p_two_orb, (mo_tot_num,mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)]
 &BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)]
 &BEGIN_PROVIDER [ double precision, total_corr_e_1h2p]
@@ -252,6 +308,8 @@ END_PROVIDER
      total_corr_e_1h2p += contrib
      corr_energy_1h2p_per_orb_ab(i_hole) += contrib
      corr_energy_1h2p_per_orb_ab(j_hole) += contrib
+     corr_energy_1h2p_two_orb(k_part,l_part) += contrib
+     corr_energy_1h2p_two_orb(l_part,k_part) += contrib
     enddo 
    enddo
   enddo
@@ -282,6 +340,8 @@ END_PROVIDER
      total_corr_e_1h2p += contrib
      corr_energy_1h2p_per_orb_aa(i_hole) += contrib
      corr_energy_1h2p_per_orb_ab(j_hole) += contrib
+     corr_energy_1h2p_two_orb(k_part,l_part) += contrib
+     corr_energy_1h2p_two_orb(l_part,k_part) += contrib
     enddo 
    enddo
   enddo
@@ -312,6 +372,8 @@ END_PROVIDER
      total_corr_e_1h2p += contrib
      corr_energy_1h2p_per_orb_bb(i_hole) += contrib
      corr_energy_1h2p_per_orb_ab(j_hole) += contrib
+     corr_energy_1h2p_two_orb(k_part,l_part) += contrib
+     corr_energy_1h2p_two_orb(l_part,k_part) += contrib
     enddo 
    enddo
   enddo

plugins/FOBOCI/dress_simple.irp.f

@@ -58,24 +58,24 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen
        call i_h_j(det_buffer(1,1,i),det_buffer(1,1,i),Nint,haa)
        f = 1.d0/(E_ref-haa)
        
-     if(second_order_h)then
+!    if(second_order_h)then
       lambda_i = f
-     else
+!    else
-       ! You write the new Hamiltonian matrix
+!      ! You write the new Hamiltonian matrix
-       do k = 1, Ndet_generators
+!      do k = 1, Ndet_generators
-        H_matrix_tmp(k,Ndet_generators+1) = H_array(k)
+!       H_matrix_tmp(k,Ndet_generators+1) = H_array(k)
-        H_matrix_tmp(Ndet_generators+1,k) = H_array(k)
+!       H_matrix_tmp(Ndet_generators+1,k) = H_array(k)
-       enddo
+!      enddo
-       H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa
+!      H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa
-       ! Then diagonalize it
+!      ! Then diagonalize it
-       call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1)
+!      call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1)
-       ! Then you extract the effective denominator
+!      ! Then you extract the effective denominator
-       accu = 0.d0
+!      accu = 0.d0
-       do k = 1, Ndet_generators
+!      do k = 1, Ndet_generators
-        accu += eigenvectors(k,1) * H_array(k)
+!       accu += eigenvectors(k,1) * H_array(k)
-       enddo
+!      enddo
-       lambda_i = eigenvectors(Ndet_generators+1,1)/accu
+!      lambda_i = eigenvectors(Ndet_generators+1,1)/accu
-     endif
+!    endif
      do k=1,idx(0)
        contrib = H_array(idx(k)) * H_array(idx(k)) * lambda_i
        delta_ij_generators_(idx(k), idx(k)) += contrib

plugins/FOBOCI/fobo_scf.irp.f

@@ -1,8 +1,8 @@
 program foboscf
  implicit none
+ call run_prepare
  no_oa_or_av_opt = .True.
  touch no_oa_or_av_opt
- call run_prepare
  call routine_fobo_scf
  call save_mos
 
@@ -10,6 +10,8 @@ end
 
 subroutine run_prepare
  implicit none
+  no_oa_or_av_opt = .False.
+  touch no_oa_or_av_opt
   call damping_SCF
   call diag_inactive_virt_and_update_mos
 end
@@ -22,6 +24,28 @@ subroutine routine_fobo_scf
  character*(64) :: label
  label = "Natural"
  do i = 1, 5
+  print*,'*******************************************************************************'
+  print*,'*******************************************************************************'
+  print*,'FOBO-SCF Iteration ',i
+  print*,'*******************************************************************************'
+  print*,'*******************************************************************************'
+  if(speed_up_convergence_foboscf)then
+   if(i==3)then
+    threshold_lmct = max(threshold_lmct,0.001)
+    threshold_mlct = max(threshold_mlct,0.05)
+    soft_touch threshold_lmct threshold_mlct
+   endif
+   if(i==4)then
+    threshold_lmct = max(threshold_lmct,0.005)
+    threshold_mlct = max(threshold_mlct,0.07)
+    soft_touch threshold_lmct threshold_mlct
+   endif
+   if(i==5)then
+    threshold_lmct = max(threshold_lmct,0.01)
+    threshold_mlct = max(threshold_mlct,0.1)
+    soft_touch threshold_lmct threshold_mlct
+   endif
+  endif
   call FOBOCI_lmct_mlct_old_thr
   call save_osoci_natural_mos
   call damping_SCF

plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f

@@ -9,12 +9,9 @@ subroutine FOBOCI_lmct_mlct_old_thr
  double precision :: norm_tmp(N_states),norm_total(N_states)
  logical :: test_sym
  double precision :: thr,hij
- double precision :: threshold
  double precision, allocatable :: dressing_matrix(:,:)
  logical :: verbose,is_ok
  verbose = .True.
- threshold = threshold_singles
- print*,'threshold = ',threshold
  thr = 1.d-12
  allocate(unpaired_bitmask(N_int,2))
  allocate (occ(N_int*bit_kind_size,2))
@@ -36,11 +33,14 @@ subroutine FOBOCI_lmct_mlct_old_thr
  print*,''
  print*,''
  print*,'DOING FIRST LMCT !!'
+ print*,'Threshold_lmct = ',threshold_lmct
  integer(bit_kind) , allocatable :: zero_bitmask(:,:)
  integer(bit_kind) , allocatable :: psi_singles(:,:,:)
+ logical :: lmct
  double precision, allocatable :: psi_singles_coef(:,:)
  allocate( zero_bitmask(N_int,2) )
   do i = 1, n_inact_orb
+   lmct = .True.
    integer :: i_hole_osoci
    i_hole_osoci = list_inact(i)
    print*,'--------------------------'
@@ -55,7 +55,7 @@ subroutine FOBOCI_lmct_mlct_old_thr
    print*,'Passed set generators'
    call set_bitmask_particl_as_input(reunion_of_bitmask)
    call set_bitmask_hole_as_input(reunion_of_bitmask)
-   call is_a_good_candidate(threshold,is_ok,verbose)
+   call is_a_good_candidate(threshold_lmct,is_ok,verbose)
    print*,'is_ok = ',is_ok
    if(.not.is_ok)cycle
    allocate(dressing_matrix(N_det_generators,N_det_generators))
@@ -81,6 +81,9 @@ subroutine FOBOCI_lmct_mlct_old_thr
     call set_bitmask_particl_as_input(reunion_of_bitmask)
     call set_bitmask_hole_as_input(reunion_of_bitmask)
     call all_single
+!   if(dressing_2h2p)then
+!    call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct)
+!   endif
 
 !   ! Change the mask of the holes and particles to perform all the 
 !   ! double excitations that starts from the active space in order 
@@ -148,9 +151,12 @@ subroutine FOBOCI_lmct_mlct_old_thr
  if(.True.)then
   print*,''
   print*,'DOING THEN THE MLCT !!'
+  print*,'Threshold_mlct = ',threshold_mlct
+   lmct = .False.
    do i = 1, n_virt_orb
     integer :: i_particl_osoci
     i_particl_osoci = list_virt(i)
+
     print*,'--------------------------'
     ! First set the current generators to the one of restart
     call set_generators_to_generators_restart
@@ -172,7 +178,7 @@ subroutine FOBOCI_lmct_mlct_old_thr
     call set_bitmask_particl_as_input(reunion_of_bitmask)
     call set_bitmask_hole_as_input(reunion_of_bitmask)
 !!  ! so all the mono excitation on the new generators 
-    call is_a_good_candidate(threshold,is_ok,verbose)
+    call is_a_good_candidate(threshold_mlct,is_ok,verbose)
     print*,'is_ok = ',is_ok
     if(.not.is_ok)cycle
      allocate(dressing_matrix(N_det_generators,N_det_generators))
@@ -187,6 +193,9 @@ subroutine FOBOCI_lmct_mlct_old_thr
  !   call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
  !   call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
      call all_single
+!    if(dressing_2h2p)then
+!     call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct)
+!    endif
     endif
     call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
     do k = 1, N_states
@@ -221,10 +230,8 @@ subroutine FOBOCI_mlct_old
  double precision :: norm_tmp,norm_total
  logical :: test_sym
  double precision :: thr
- double precision :: threshold
  logical :: verbose,is_ok
  verbose = .False.
- threshold = 1.d-2
  thr = 1.d-12
  allocate(unpaired_bitmask(N_int,2))
  allocate (occ(N_int*bit_kind_size,2))
@@ -263,7 +270,7 @@ subroutine FOBOCI_mlct_old
    call set_bitmask_particl_as_input(reunion_of_bitmask)
    call set_bitmask_hole_as_input(reunion_of_bitmask)
 !  ! so all the mono excitation on the new generators 
-   call is_a_good_candidate(threshold,is_ok,verbose)
+   call is_a_good_candidate(threshold_mlct,is_ok,verbose)
    print*,'is_ok = ',is_ok
    is_ok =.True.
    if(.not.is_ok)cycle
@@ -297,10 +304,8 @@ subroutine FOBOCI_lmct_old
  double precision :: norm_tmp,norm_total
  logical :: test_sym
  double precision :: thr
- double precision :: threshold
  logical :: verbose,is_ok
  verbose = .False.
- threshold = 1.d-2
  thr = 1.d-12
  allocate(unpaired_bitmask(N_int,2))
  allocate (occ(N_int*bit_kind_size,2))
@@ -337,7 +342,7 @@ subroutine FOBOCI_lmct_old
    call set_generators_to_psi_det
    call set_bitmask_particl_as_input(reunion_of_bitmask)
    call set_bitmask_hole_as_input(reunion_of_bitmask)
-   call is_a_good_candidate(threshold,is_ok,verbose)
+   call is_a_good_candidate(threshold_lmct,is_ok,verbose)
    print*,'is_ok = ',is_ok
    if(.not.is_ok)cycle
 !  ! so all the mono excitation on the new generators 

plugins/FOBOCI/new_approach.irp.f

@@ -46,7 +46,7 @@ subroutine new_approach
 
 
  verbose = .True.
- threshold = threshold_singles
+ threshold = threshold_lmct
  print*,'threshold = ',threshold
  thr = 1.d-12
  print*,''
@@ -623,14 +623,14 @@ subroutine new_approach
   print*,'ACTIVE ORBITAL  ',iorb
   do j = 1, n_inact_orb
    jorb = list_inact(j)
-   if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_lmct)then
     print*,'INACTIVE  '
     print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb))
    endif
   enddo
   do j = 1, n_virt_orb
    jorb = list_virt(j)
-   if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_mlct)then
     print*,'VIRT      '
     print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb))
    endif

plugins/FOBOCI/routines_foboci.irp.f

@@ -387,14 +387,14 @@ subroutine save_osoci_natural_mos
   print*,'ACTIVE ORBITAL  ',iorb
   do j = 1, n_inact_orb
    jorb = list_inact(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then
     print*,'INACTIVE  '
     print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
   enddo
   do j = 1, n_virt_orb
    jorb = list_virt(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then
     print*,'VIRT      '
     print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
@@ -519,14 +519,14 @@ subroutine set_osoci_natural_mos
   print*,'ACTIVE ORBITAL  ',iorb
   do j = 1, n_inact_orb
    jorb = list_inact(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then
     print*,'INACTIVE  '
     print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
   enddo
   do j = 1, n_virt_orb
    jorb = list_virt(j)
-   if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
+   if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then
     print*,'VIRT      '
     print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
@@ -607,3 +607,206 @@ end
    call print_hcc
  end
 
+
+
+ subroutine dress_diag_elem_2h1p(dressing_H_mat_elem,ndet,lmct,i_hole)
+  use bitmasks
+  double precision, intent(inout) :: dressing_H_mat_elem(Ndet)
+  integer, intent(in) :: ndet,i_hole
+  logical, intent(in) :: lmct
+  ! if lmct = .True. ===> LMCT
+  ! else             ===> MLCT
+  implicit none
+  integer :: i
+  integer :: n_p,n_h,number_of_holes,number_of_particles
+  integer :: exc(0:2,2,2)
+  integer :: degree
+  double precision :: phase
+  integer :: h1,h2,p1,p2,s1,s2
+  do i = 1, N_det
+
+   n_h = number_of_holes(psi_det(1,1,i))
+   n_p = number_of_particles(psi_det(1,1,i))
+   call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int)
+   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+   if (n_h == 0.and.n_p==0)then ! CAS
+    dressing_H_mat_elem(i)+= total_corr_e_2h1p
+    if(lmct)then
+     dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(i_hole) - corr_energy_2h1p_per_orb_bb(i_hole)    
+    endif
+   endif
+   if (n_h == 1.and.n_p==0)then ! 1h
+     dressing_H_mat_elem(i)+= 0.d0
+   else if (n_h == 0.and.n_p==1)then ! 1p
+    dressing_H_mat_elem(i)+= total_corr_e_2h1p
+    dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(p1) - corr_energy_2h1p_per_orb_aa(p1)  
+   else if (n_h == 1.and.n_p==1)then ! 1h1p
+!   if(degree==1)then
+     dressing_H_mat_elem(i)+= total_corr_e_2h1p
+     dressing_H_mat_elem(i)+= - corr_energy_2h1p_per_orb_ab(h1) 
+!   else 
+!    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+!    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+!    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) &  
+!    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2))
+!    dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1))
+!   endif
+   else if (n_h == 2.and.n_p==1)then ! 2h1p
+    dressing_H_mat_elem(i)+= 0.d0
+   else if (n_h == 1.and.n_p==2)then ! 1h2p
+    dressing_H_mat_elem(i)+= total_corr_e_2h1p
+    dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(h1) 
+   endif
+  enddo
+ 
+ end
+
+ subroutine dress_diag_elem_1h2p(dressing_H_mat_elem,ndet,lmct,i_hole)
+  use bitmasks
+  double precision, intent(inout) :: dressing_H_mat_elem(Ndet)
+  integer, intent(in) :: ndet,i_hole
+  logical, intent(in) :: lmct
+  ! if lmct = .True. ===> LMCT
+  ! else             ===> MLCT
+  implicit none
+  integer :: i
+  integer :: n_p,n_h,number_of_holes,number_of_particles
+  integer :: exc(0:2,2,2)
+  integer :: degree
+  double precision :: phase
+  integer :: h1,h2,p1,p2,s1,s2
+  do i = 1, N_det
+
+   n_h = number_of_holes(psi_det(1,1,i))
+   n_p = number_of_particles(psi_det(1,1,i))
+   call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int)
+   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+   if (n_h == 0.and.n_p==0)then ! CAS
+    dressing_H_mat_elem(i)+= total_corr_e_1h2p
+    if(.not.lmct)then
+     dressing_H_mat_elem(i) += - corr_energy_1h2p_per_orb_ab(i_hole) - corr_energy_1h2p_per_orb_aa(i_hole)    
+    endif
+   endif
+   if (n_h == 1.and.n_p==0)then ! 1h
+     dressing_H_mat_elem(i)+= total_corr_e_1h2p - corr_energy_1h2p_per_orb_ab(h1)
+   else if (n_h == 0.and.n_p==1)then ! 1p
+    dressing_H_mat_elem(i)+= 0.d0
+   else if (n_h == 1.and.n_p==1)then ! 1h1p
+    if(degree==1)then
+     dressing_H_mat_elem(i)+= total_corr_e_1h2p
+     dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) 
+    else 
+     dressing_H_mat_elem(i) +=0.d0
+    endif
+!    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+!    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+!    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) &  
+!    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2))
+!    dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1))
+!   endif
+   else if (n_h == 2.and.n_p==1)then ! 2h1p
+    dressing_H_mat_elem(i)+= total_corr_e_1h2p
+    dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) - corr_energy_1h2p_per_orb_ab(h1)
+   else if (n_h == 1.and.n_p==2)then ! 1h2p
+    dressing_H_mat_elem(i) += 0.d0
+   endif
+  enddo
+ 
+ end
+
+ subroutine dress_diag_elem_2h2p(dressing_H_mat_elem,ndet)
+  use bitmasks
+  double precision, intent(inout) :: dressing_H_mat_elem(Ndet)
+  integer, intent(in) :: ndet
+  implicit none
+  integer :: i
+  integer :: n_p,n_h,number_of_holes,number_of_particles
+  integer :: exc(0:2,2,2)
+  integer :: degree
+  double precision :: phase
+  integer :: h1,h2,p1,p2,s1,s2
+  do i = 1, N_det
+   dressing_H_mat_elem(i)+= total_corr_e_2h2p
+
+   n_h = number_of_holes(psi_det(1,1,i))
+   n_p = number_of_particles(psi_det(1,1,i))
+   call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int)
+   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+   if (n_h == 1.and.n_p==0)then ! 1h
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+   else if (n_h == 0.and.n_p==1)then ! 1p
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) &  
+    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1))
+   else if (n_h == 1.and.n_p==1)then ! 1h1p
+    if(degree==1)then
+     dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+     - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+     dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) &  
+     - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1))
+     dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_a(h1,p1) + corr_energy_2h2p_for_1h1p_b(h1,p1))
+    else 
+     dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+     - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+     dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) &  
+     - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2))
+     dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1))
+    endif
+   else if (n_h == 2.and.n_p==1)then ! 2h1p
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1) & 
+                              - corr_energy_2h2p_per_orb_ab(h2)                                   & 
+                              - 0.5d0 * ( corr_energy_2h2p_per_orb_bb(h2) + corr_energy_2h2p_per_orb_bb(h2)) 
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1)
+    if(s1.ne.s2)then
+      dressing_H_mat_elem(i) +=  corr_energy_2h2p_ab_2_orb(h1,h2)
+    else 
+      dressing_H_mat_elem(i) +=  corr_energy_2h2p_bb_2_orb(h1,h2)
+    endif
+   else if (n_h == 1.and.n_p==2)then ! 1h2p
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) &  
+    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1))
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) &  
+    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1))
+    dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) &  
+    - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2))
+    if(s1.ne.s2)then
+      dressing_H_mat_elem(i) +=  corr_energy_2h2p_ab_2_orb(p1,p2)
+    else 
+      dressing_H_mat_elem(i) +=  corr_energy_2h2p_bb_2_orb(p1,p2)
+    endif
+   endif
+  enddo
+ 
+ end
+
+ subroutine diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole,lmct)
+ implicit none
+ double precision, allocatable :: dressing_H_mat_elem(:),energies(:)
+  integer, intent(in) :: i_hole
+  logical, intent(in) :: lmct
+  ! if lmct = .True. ===> LMCT
+  ! else             ===> MLCT
+ integer :: i
+ double precision :: hij
+ allocate(dressing_H_mat_elem(N_det),energies(N_states_diag))
+  print*,''
+  print*,'dressing with the 2h2p in a CC logic'
+  print*,''
+  do i = 1, N_det
+   call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
+   dressing_H_mat_elem(i) = hij
+  enddo
+  call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det)
+  call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,i_hole)
+  call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,i_hole)
+  call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states_diag,N_int,output_determinants)
+  do i = 1, 2
+   print*,'psi_coef = ',psi_coef(i,1)
+  enddo
+  
+  
+ deallocate(dressing_H_mat_elem)
+
+
+
+ end