Merge pull request #156 from eginer/master

FOBO-SCF algorithm in a stable version,

config/ifort.cfg

@@ -51,7 +51,7 @@ FCFLAGS  : -xSSE4.2 -O2 -ip -ftz
 #
 [DEBUG]
 FC      : -g -traceback
-FCFLAGS : -xSSE2 -C  -fpe0 
+FCFLAGS : -xSSE2 -C  
 IRPF90_FLAGS : --openmp
 
 # OpenMP flags

plugins/All_singles/H_apply.irp.f

@@ -8,10 +8,9 @@ s.unset_skip()
 s.filter_only_1h1p()
 print s
 
-s = H_apply("just_mono")
+s = H_apply("just_mono",do_double_exc=False)
 s.set_selection_pt2("epstein_nesbet_2x2")
 s.unset_skip()
-s.unset_double_excitations()
 print s
 
 END_SHELL

plugins/All_singles/all_singles.irp.f

@@ -15,7 +15,7 @@ subroutine routine
   integer                        :: N_st, degree
   double precision,allocatable :: E_before(:)
   integer :: n_det_before
-  N_st = N_states
+  N_st = N_states_diag
   allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
   i = 0
   print*,'N_det = ',N_det

plugins/CAS_SD/H_apply.irp.f

@@ -20,22 +20,18 @@ print s
 
 
 s = H_apply("CAS_S",do_double_exc=False)
-s.unset_double_excitations()
 print s
 
 s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
-s.unset_double_excitations()
 s.set_selection_pt2("epstein_nesbet_2x2")
 s.unset_skip()
 print s
 
 s = H_apply("CAS_S_selected",do_double_exc=False)
-s.unset_double_excitations()
 s.set_selection_pt2("epstein_nesbet_2x2")
 print s
 
 s = H_apply("CAS_S_PT2",do_double_exc=False)
-s.unset_double_excitations()
 s.set_perturbation("epstein_nesbet_2x2")
 print s
 

plugins/DDCI_selected/ddci.irp.f

@@ -3,10 +3,10 @@ program ddci
   integer                        :: i,k
 
   
-  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:)
   integer                        :: N_st, degree
-  N_st = N_states
+  N_st = N_states_diag
-  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
   character*(64)                 :: perturbation
   
   pt2 = 1.d0
@@ -27,6 +27,8 @@ program ddci
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
   endif
+   call set_bitmask_particl_as_input(reunion_of_bitmask)
+   call set_bitmask_hole_as_input(reunion_of_bitmask)
 
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
     call H_apply_DDCI_selection(pt2, norm_pert, H_pert_diag,  N_st)
@@ -47,8 +49,21 @@ program ddci
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
     print *,  'E        = ', CI_energy
-    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  'E+PT2    = ', E_before+pt2
     print *,  '-----'
+    if(N_states_diag.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
+     enddo
+    endif
+    E_before = CI_energy
     call ezfio_set_ddci_selected_energy(CI_energy)
   enddo
   if(do_pt2_end)then

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/H_apply.irp.f

@@ -18,8 +18,22 @@ print s
 
 
 
-s = H_apply("standard")
+s = H_apply("only_1h2p")
 s.set_selection_pt2("epstein_nesbet")
+s.filter_only_1h2p()
+s.unset_skip()
+print s
+
+s = H_apply("only_2h2p")
+s.set_selection_pt2("epstein_nesbet")
+s.filter_only_2h2p()
+s.unset_skip()
+print s
+
+
+s = H_apply("only_2p")
+s.set_selection_pt2("epstein_nesbet")
+s.filter_only_2p()
 s.unset_skip()
 print s
 

plugins/FOBOCI/H_apply_dressed_autonom.irp.f

@@ -437,8 +437,8 @@ subroutine H_apply_dressed_pert(delta_ij_generators_,  Ndet_generators,psi_det_g
   
   
     integer, intent(in) :: Ndet_generators
-    integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators),E_ref
+    integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
-    double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
+    double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
 
   
   integer                        :: i_generator, nmax

plugins/FOBOCI/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Generators_restart Selectors_no_sorted
+Perturbation Selectors_no_sorted Hartree_Fock

plugins/FOBOCI/all_singles.irp.f

@@ -6,9 +6,9 @@ subroutine all_single
   double precision,allocatable :: E_before(:)
   N_st = N_states
   allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
-  selection_criterion = 1.d-8
+  selection_criterion = 0.d0
   soft_touch selection_criterion
-  threshold_davidson = 1.d-5
+  threshold_davidson = 1.d-9
   soft_touch threshold_davidson davidson_criterion
   i = 0
   print*,'Doing all the mono excitations !'
@@ -52,10 +52,173 @@ subroutine all_single
      enddo
     endif
     E_before = CI_energy
+    !!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO
+    exit
   enddo
-  threshold_davidson = 1.d-10
+! threshold_davidson = 1.d-8
+! soft_touch threshold_davidson davidson_criterion
+! call diagonalize_CI
+  print*,'Final Step '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+  do i = 1, max(2,N_det_generators)
+   print*,'psi_coef = ',psi_coef(i,1)
+  enddo
+  deallocate(pt2,norm_pert,E_before)
+end
+
+subroutine all_1h2p
+  implicit none
+  integer                        :: i,k
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  double precision,allocatable :: E_before(:)
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
+  selection_criterion = 0.d0
+  soft_touch selection_criterion
+  threshold_davidson = 1.d-5
   soft_touch threshold_davidson davidson_criterion
-  call diagonalize_CI
+  i = 0
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,'*****************************'
+  print*,'Doing all the 1h2P excitations'
+  print*,'*****************************'
+  print*,''
+  print*,''
+  print*,'N_det = ',N_det
+  print*,'n_det_max = ',n_det_max
+  print*,'pt2_max = ',pt2_max
+  print*,'N_det_generators = ',N_det_generators
+  pt2=-1.d0
+  E_before = ref_bitmask_energy
+ 
+  print*,'Initial Step '
+  print*,'Inital determinants '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+  n_det_max = 100000
+  i = 0
+  do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    i += 1
+    print*,'-----------------------'
+    print*,'i = ',i
+    call H_apply_only_1h2p(pt2, norm_pert, H_pert_diag,  N_st)
+    call diagonalize_CI
+    print*,'N_det = ',N_det
+    print*,'E        = ',CI_energy(1)
+    print*,'pt2      = ',pt2(1)
+    print*,'E+PT2    = ',E_before + pt2(1)
+    if(N_states_diag.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
+     enddo
+    endif
+    E_before = CI_energy
+    
+  enddo
+  print*,'Final Step '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+
+  do i = 1, 2
+   print*,'psi_coef = ',psi_coef(i,1)
+  enddo
+  deallocate(pt2,norm_pert,E_before)
+end
+
+subroutine all_2h2p
+  implicit none
+  integer                        :: i,k
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  double precision,allocatable :: E_before(:)
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
+  selection_criterion = 0.d0
+  soft_touch selection_criterion
+  threshold_davidson = 1.d-5
+  soft_touch threshold_davidson davidson_criterion
+  i = 0
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,'*****************************'
+  print*,'Doing all the 2h2P excitations'
+  print*,'*****************************'
+  print*,''
+  print*,''
+  print*,'N_det = ',N_det
+  print*,'n_det_max = ',n_det_max
+  print*,'pt2_max = ',pt2_max
+  print*,'N_det_generators = ',N_det_generators
+  pt2=-1.d0
+  E_before = ref_bitmask_energy
+ 
+  print*,'Initial Step '
+  print*,'Inital determinants '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+  n_det_max = 100000
+  i = 0
+  do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    i += 1
+    print*,'-----------------------'
+    print*,'i = ',i
+    call H_apply_only_2h2p(pt2, norm_pert, H_pert_diag,  N_st)
+    call diagonalize_CI
+    print*,'N_det = ',N_det
+    print*,'E        = ',CI_energy(1)
+    print*,'pt2      = ',pt2(1)
+    print*,'E+PT2    = ',E_before + pt2(1)
+    if(N_states_diag.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
+     enddo
+    endif
+    E_before = CI_energy
+    
+  enddo
   print*,'Final Step '
   print*,'N_det = ',N_det
   do i = 1, N_states_diag
@@ -67,10 +230,89 @@ subroutine all_single
   do i = 1, 2
    print*,'psi_coef = ',psi_coef(i,1)
   enddo
-! call save_wavefunction
   deallocate(pt2,norm_pert,E_before)
 end
 
+subroutine all_2p
+  implicit none
+  integer                        :: i,k
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  double precision,allocatable :: E_before(:)
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
+  selection_criterion = 0.d0
+  soft_touch selection_criterion
+  threshold_davidson = 1.d-5
+  soft_touch threshold_davidson davidson_criterion
+  i = 0
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,''
+  print*,'*****************************'
+  print*,'Doing all the 2P excitations'
+  print*,'*****************************'
+  print*,''
+  print*,''
+  print*,'N_det = ',N_det
+  print*,'n_det_max = ',n_det_max
+  print*,'pt2_max = ',pt2_max
+  print*,'N_det_generators = ',N_det_generators
+  pt2=-1.d0
+  E_before = ref_bitmask_energy
+ 
+  print*,'Initial Step '
+  print*,'Inital determinants '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+  n_det_max = 100000
+  i = 0
+  do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    i += 1
+    print*,'-----------------------'
+    print*,'i = ',i
+    call H_apply_only_2p(pt2, norm_pert, H_pert_diag,  N_st)
+    call diagonalize_CI
+    print*,'N_det = ',N_det
+    print*,'E        = ',CI_energy(1)
+    print*,'pt2      = ',pt2(1)
+    print*,'E+PT2    = ',E_before + pt2(1)
+    if(N_states_diag.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_st
+      print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
+     enddo
+    endif
+    E_before = CI_energy
+    
+  enddo
+  print*,'Final Step '
+  print*,'N_det = ',N_det
+  do i = 1, N_states_diag
+   print*,''
+   print*,'i = ',i
+   print*,'E        = ',CI_energy(i)
+   print*,'S^2      = ',CI_eigenvectors_s2(i)
+  enddo
+  deallocate(pt2,norm_pert,E_before)
+  do i = 1, 2
+   print*,'psi_coef = ',psi_coef(i,1)
+  enddo
+end
+
 subroutine all_single_no_1h_or_1p
   implicit none
   integer                        :: i,k
@@ -79,6 +321,8 @@ subroutine all_single_no_1h_or_1p
   double precision,allocatable :: E_before(:)
   N_st = N_states
   allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
+  selection_criterion = 0.d0
+  soft_touch selection_criterion
   threshold_davidson = 1.d-5
   soft_touch threshold_davidson davidson_criterion
   i = 0
@@ -124,7 +368,7 @@ subroutine all_single_no_1h_or_1p
     endif
     E_before = CI_energy
   enddo
-  threshold_davidson = 1.d-10
+  threshold_davidson = 1.d-16
   soft_touch threshold_davidson davidson_criterion
   call diagonalize_CI
   print*,'Final Step '
@@ -215,85 +459,6 @@ subroutine all_single_no_1h_or_1p_or_2p
   deallocate(pt2,norm_pert,E_before)
 end
 
-
-subroutine all_2p
-  implicit none
-  integer                        :: i,k
-  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
-  integer                        :: N_st, degree
-  double precision,allocatable :: E_before(:)
-  N_st = N_states
-  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
-  selection_criterion = 0.d0
-  soft_touch selection_criterion
-  threshold_davidson = 1.d-5
-  soft_touch threshold_davidson davidson_criterion
-  i = 0
-  print*,''
-  print*,''
-  print*,''
-  print*,''
-  print*,''
-  print*,'*****************************'
-  print*,'Doing all the 2P excitations'
-  print*,'*****************************'
-  print*,''
-  print*,''
-  print*,'N_det = ',N_det
-  print*,'n_det_max = ',n_det_max
-  print*,'pt2_max = ',pt2_max
-  print*,'N_det_generators = ',N_det_generators
-  pt2=-1.d0
-  E_before = ref_bitmask_energy
- 
-  print*,'Initial Step '
-  print*,'Inital determinants '
-  print*,'N_det = ',N_det
-  do i = 1, N_states_diag
-   print*,''
-   print*,'i = ',i
-   print*,'E        = ',CI_energy(i)
-   print*,'S^2      = ',CI_eigenvectors_s2(i)
-  enddo
-  n_det_max = 100000
-  i = 0
-  do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
-    i += 1
-    print*,'-----------------------'
-    print*,'i = ',i
-    call H_apply_standard(pt2, norm_pert, H_pert_diag,  N_st)
-    call diagonalize_CI
-    print*,'N_det = ',N_det
-    print*,'E        = ',CI_energy(1)
-    print*,'pt2      = ',pt2(1)
-    print*,'E+PT2    = ',E_before + pt2(1)
-    if(N_states_diag.gt.1)then
-     print*,'Variational Energy difference'
-     do i = 2, N_st
-      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
-     enddo
-    endif
-    if(N_states.gt.1)then
-     print*,'Variational + perturbative Energy difference'
-     do i = 2, N_st
-      print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
-     enddo
-    endif
-    E_before = CI_energy
-    
-  enddo
-  print*,'Final Step '
-  print*,'N_det = ',N_det
-  do i = 1, N_states_diag
-   print*,''
-   print*,'i = ',i
-   print*,'E        = ',CI_energy(i)
-   print*,'S^2      = ',CI_eigenvectors_s2(i)
-  enddo
-! call save_wavefunction
-  deallocate(pt2,norm_pert,E_before)
-end
-
 subroutine all_1h_1p_routine
   implicit none
   integer                        :: i,k

plugins/FOBOCI/all_singles_split.irp.f

@@ -5,7 +5,7 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N
  integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators_input)
  double precision, intent(inout) :: dressing_matrix(Ndet_generators_input,Ndet_generators_input)
  double precision, intent(in)    :: psi_coef_generators_input(ndet_generators_input,n_states)
- integer :: i,i_hole
+ integer :: i,i_hole,j
  n_det_max_jacobi = 50
  soft_touch n_det_max_jacobi
  do i = 1, n_inact_orb
@@ -22,56 +22,339 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N
    call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
    call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
    call all_single
-   threshold_davidson = 1.d-10
+!  call diagonalize_CI_SC2
-   soft_touch threshold_davidson davidson_criterion
+!  call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
-   call diagonalize_CI
    call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
  enddo
+
+ do i = 1, n_act_orb
+   i_hole = list_act(i)
+   print*,''
+   print*,'Doing all the single excitations from the orbital '
+   print*,i_hole
+   print*,''
+   print*,''
+   threshold_davidson = 1.d-4
+   soft_touch threshold_davidson davidson_criterion
+   call modify_bitmasks_for_hole(i_hole)
+   call set_bitmask_particl_as_input(reunion_of_bitmask)
+   call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
+   call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
+   call all_single
+!  call diagonalize_CI_SC2
+!  call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
+   call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
+ enddo
+
+ do i = 1, n_virt_orb
+   i_hole = list_virt(i)
+   print*,''
+   print*,'Doing all the single excitations from the orbital '
+   print*,i_hole
+   print*,''
+   print*,''
+   threshold_davidson = 1.d-4
+   soft_touch threshold_davidson davidson_criterion
+   call modify_bitmasks_for_hole(i_hole)
+   call set_bitmask_particl_as_input(reunion_of_bitmask)
+   call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
+   call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
+   call all_single
+!  call diagonalize_CI_SC2
+!  call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
+   call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
+ enddo
+
  n_det_max_jacobi = 1000
  soft_touch n_det_max_jacobi
 end
 
 
-subroutine all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
+
+subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
  implicit none
  use bitmasks
+ integer, intent(in)             :: i_particl
  double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
- double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators)
+ double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators)
- integer :: i,i_hole
+ double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
+ integer :: i,j
  n_det_max_jacobi = 50
  soft_touch n_det_max_jacobi
 
- integer :: n_det_1h1p,n_det_2h1p
+ call all_single
- integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
- integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
- integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
- double precision, allocatable :: psi_ref_coef_out(:,:)
- double precision, allocatable :: psi_coef_1h1p(:,:)
- double precision, allocatable :: psi_coef_2h1p(:,:)
- call all_single_no_1h_or_1p
  
  threshold_davidson = 1.d-12
  soft_touch threshold_davidson davidson_criterion
  call diagonalize_CI
- call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p)
- allocate(psi_ref_out(N_int,2,N_det_generators))
- allocate(psi_1h1p(N_int,2,n_det_1h1p))
- allocate(psi_2h1p(N_int,2,n_det_2h1p))
- allocate(psi_ref_coef_out(N_det_generators,N_states))
- allocate(psi_coef_1h1p(n_det_1h1p,N_states))
- allocate(psi_coef_2h1p(n_det_2h1p,N_states))
- call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p)
- call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
-                                                        psi_1h1p,psi_coef_1h1p,n_det_1h1p)
- call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
-                                                        psi_2h1p,psi_coef_2h1p,n_det_2h1p)
 
- deallocate(psi_ref_out)
+
- deallocate(psi_1h1p)
+
- deallocate(psi_2h1p)
+ double precision, allocatable  :: matrix_ref_1h_1p(:,:)
- deallocate(psi_ref_coef_out)
+ double precision, allocatable  :: matrix_ref_1h_1p_dressing_1h1p(:,:)
- deallocate(psi_coef_1h1p)
+ double precision, allocatable  :: matrix_ref_1h_1p_dressing_1h2p(:,:)
- deallocate(psi_coef_2h1p)
+ double precision, allocatable  :: psi_coef_ref_1h_1p(:,:)
+ double precision, allocatable  :: psi_coef_1h1p(:,:)
+ double precision, allocatable  :: psi_coef_1h2p(:,:)
+ integer(bit_kind), allocatable :: psi_det_1h2p(:,:,:)
+ integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:)
+ integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:)
+ integer :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
+ double precision :: hka
+ double precision,allocatable :: eigenvectors(:,:), eigenvalues(:)
+
+
+ call give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p)
+
+ allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(matrix_ref_1h_1p_dressing_1h2p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states))
+ allocate(psi_det_1h2p(N_int,2,n_det_1h2p), psi_coef_1h2p(n_det_1h2p,N_states))
+ allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states))
+
+ call give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
+                                                          psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p)
+
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, n_det_ref_1h_1p
+   call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka)
+   matrix_ref_1h_1p(i,j) = hka
+  enddo
+ enddo
+ matrix_ref_1h_1p_dressing_1h1p = 0.d0
+ matrix_ref_1h_1p_dressing_1h2p = 0.d0
+ call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h2p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
+                                                        psi_det_1h2p,psi_coef_1h2p,n_det_1h2p)
+ call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
+                                                        psi_det_1h1p,psi_coef_1h1p,n_det_1h1p)
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, n_det_ref_1h_1p
+   matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_1h2p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j)
+  enddo
+ enddo
+ 
+ allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p))
+ call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p)
+!do j = 1, n_det_ref_1h_1p
+! print*,'coef = ',eigenvectors(j,1)
+!enddo
+ print*,''
+ print*,'-----------------------'
+ print*,'-----------------------'
+ print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion
+ print*,'-----------------------'
+ ! Extract the 
+ integer, allocatable :: index_generator(:)
+ integer :: n_det_generators_tmp,degree
+ n_det_generators_tmp = 0
+ allocate(index_generator(n_det_ref_1h_1p))
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, N_det_generators
+   call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(1,1,i), degree, N_int)
+   if(degree == 0)then
+    n_det_generators_tmp +=1
+    index_generator(n_det_generators_tmp) = i
+   endif
+  enddo
+ enddo
+ if(n_det_generators_tmp .ne. n_det_generators)then
+  print*,'PB !!!'
+  print*,'if(n_det_generators_tmp .ne. n_det_genrators)then'
+  stop
+ endif
+ do i = 1, N_det_generators
+  print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1)
+  do j = 1, N_det_generators
+   dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j))
+   dressing_matrix_1h2p(i,j) += matrix_ref_1h_1p_dressing_1h2p(index_generator(i),index_generator(j))
+  enddo
+ enddo
+ print*,'-----------------------'
+ print*,'-----------------------'
+ 
+
+ deallocate(matrix_ref_1h_1p)
+ deallocate(matrix_ref_1h_1p_dressing_1h1p)
+ deallocate(matrix_ref_1h_1p_dressing_1h2p)
+ deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p)
+ deallocate(psi_det_1h2p, psi_coef_1h2p)
+ deallocate(psi_det_1h1p, psi_coef_1h1p)
+ deallocate(eigenvectors,eigenvalues)
+ deallocate(index_generator)
+
+
+end
+
+subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p)
+ implicit none
+ use bitmasks
+ integer, intent(in)             :: i_hole
+ double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
+ double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators)
+ double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
+ integer :: i,j
+ n_det_max_jacobi = 50
+ soft_touch n_det_max_jacobi
+
+ call all_single
+ 
+ threshold_davidson = 1.d-12
+ soft_touch threshold_davidson davidson_criterion
+ call diagonalize_CI
+
+
+
+ double precision, allocatable  :: matrix_ref_1h_1p(:,:)
+ double precision, allocatable  :: matrix_ref_1h_1p_dressing_1h1p(:,:)
+ double precision, allocatable  :: matrix_ref_1h_1p_dressing_2h1p(:,:)
+ double precision, allocatable  :: psi_coef_ref_1h_1p(:,:)
+ double precision, allocatable  :: psi_coef_1h1p(:,:)
+ double precision, allocatable  :: psi_coef_2h1p(:,:)
+ integer(bit_kind), allocatable :: psi_det_2h1p(:,:,:)
+ integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:)
+ integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:)
+ integer :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
+ double precision :: hka
+ double precision,allocatable :: eigenvectors(:,:), eigenvalues(:)
+
+
+ call give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p)
+
+ allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(matrix_ref_1h_1p_dressing_2h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
+ allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states))
+ allocate(psi_det_2h1p(N_int,2,n_det_2h1p), psi_coef_2h1p(n_det_2h1p,N_states))
+ allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states))
+
+ call give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
+                                                          psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p)
+
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, n_det_ref_1h_1p
+   call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka)
+   matrix_ref_1h_1p(i,j) = hka
+  enddo
+ enddo
+ matrix_ref_1h_1p_dressing_1h1p = 0.d0
+ matrix_ref_1h_1p_dressing_2h1p = 0.d0
+ call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_2h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
+                                                        psi_det_2h1p,psi_coef_2h1p,n_det_2h1p)
+ call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
+                                                        psi_det_1h1p,psi_coef_1h1p,n_det_1h1p)
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, n_det_ref_1h_1p
+   matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j)
+  enddo
+ enddo
+ 
+ allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p))
+ call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p)
+!do j = 1, n_det_ref_1h_1p
+! print*,'coef = ',eigenvectors(j,1)
+!enddo
+ print*,''
+ print*,'-----------------------'
+ print*,'-----------------------'
+ print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion
+ print*,'-----------------------'
+ ! Extract the 
+ integer, allocatable :: index_generator(:)
+ integer :: n_det_generators_tmp,degree
+ n_det_generators_tmp = 0
+ allocate(index_generator(n_det_ref_1h_1p))
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, N_det_generators
+   call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(1,1,i), degree, N_int)
+   if(degree == 0)then
+    n_det_generators_tmp +=1
+    index_generator(n_det_generators_tmp) = i
+   endif
+  enddo
+ enddo
+ if(n_det_generators_tmp .ne. n_det_generators)then
+  print*,'PB !!!'
+  print*,'if(n_det_generators_tmp .ne. n_det_genrators)then'
+  stop
+ endif
+ do i = 1, N_det_generators
+  print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1)
+  do j = 1, N_det_generators
+   dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j))
+   dressing_matrix_2h1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(index_generator(i),index_generator(j))
+  enddo
+ enddo
+ print*,'-----------------------'
+ print*,'-----------------------'
+ 
+
+ deallocate(matrix_ref_1h_1p)
+ deallocate(matrix_ref_1h_1p_dressing_1h1p)
+ deallocate(matrix_ref_1h_1p_dressing_2h1p)
+ deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p)
+ deallocate(psi_det_2h1p, psi_coef_2h1p)
+ deallocate(psi_det_1h1p, psi_coef_1h1p)
+ deallocate(eigenvectors,eigenvalues)
+ deallocate(index_generator)
+!return
+!
+
+!integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
+!integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
+!integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
+!integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
+!double precision, allocatable :: psi_ref_coef_out(:,:)
+!double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
+ 
+!call all_single_no_1h_or_1p
+
+!call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p)
+!allocate(psi_ref_out(N_int,2,N_det_generators))
+!allocate(psi_1h1p(N_int,2,n_det_1h1p))
+!allocate(psi_2h1p(N_int,2,n_det_2h1p))
+!allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
+!allocate(psi_ref_coef_out(N_det_generators,N_states))
+!allocate(psi_coef_1h1p(n_det_1h1p,N_states))
+!allocate(psi_coef_2h1p(n_det_2h1p,N_states))
+!allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
+!call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
+!do i = 1, n_det_extra_1h_or_1p
+!  print*,'----'
+!  print*,'c = ',psi_coef_extra_1h_or_1p(i,1)
+!  call debug_det(psi_extra_1h_or_1p(1,1,i),N_int)
+!  print*,'----'
+!enddo
+!call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
+!                                                       psi_1h1p,psi_coef_1h1p,n_det_1h1p)
+!print*,'Dressing 1h1p '
+!do j =1, N_det_generators
+! print*,' dressing ',dressing_matrix_1h1p(j,:)
+!enddo
+
+!call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
+!                                                       psi_2h1p,psi_coef_2h1p,n_det_2h1p)
+!print*,'Dressing 2h1p '
+!do j =1, N_det_generators
+! print*,' dressing ',dressing_matrix_2h1p(j,:)
+!enddo
+
+!call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
+!                                                       psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
+!print*,',dressing_matrix_extra_1h_or_1p'
+!do j =1, N_det_generators
+! print*,' dressing ',dressing_matrix_extra_1h_or_1p(j,:)
+!enddo
+
+
+!deallocate(psi_ref_out)
+!deallocate(psi_1h1p)
+!deallocate(psi_2h1p)
+!deallocate(psi_extra_1h_or_1p)
+!deallocate(psi_ref_coef_out)
+!deallocate(psi_coef_1h1p)
+!deallocate(psi_coef_2h1p)
+!deallocate(psi_coef_extra_1h_or_1p)
 
 end
 
@@ -197,47 +480,56 @@ subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
  soft_touch n_det_max_jacobi
 end
 
-subroutine all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
+!  subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
- implicit none
+!  implicit none
- use bitmasks
+!  use bitmasks
- double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
+!  integer, intent(in )            :: i_particl
- double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators)
+!  double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
- integer :: i,i_hole
+!  double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators)
- n_det_max_jacobi = 50
+!  double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
- soft_touch n_det_max_jacobi
+!  integer :: i
-
+!  n_det_max_jacobi = 50
- integer :: n_det_1h1p,n_det_1h2p
+!  soft_touch n_det_max_jacobi
- integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
+! 
- integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
+!  integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p
- integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
+!  integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
- double precision, allocatable :: psi_ref_coef_out(:,:)
+!  integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
- double precision, allocatable :: psi_coef_1h1p(:,:)
+!  integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
- double precision, allocatable :: psi_coef_1h2p(:,:)
+!  integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
- call all_single_no_1h_or_1p_or_2p
+!  double precision, allocatable :: psi_ref_coef_out(:,:)
- 
+!  double precision, allocatable :: psi_coef_1h1p(:,:)
- threshold_davidson = 1.d-12
+!  double precision, allocatable :: psi_coef_1h2p(:,:)
- soft_touch threshold_davidson davidson_criterion
+!  double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
- call diagonalize_CI
+!!!!call all_single_no_1h_or_1p_or_2p
- call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p)
+!  call all_single
- allocate(psi_ref_out(N_int,2,N_det_generators))
+!  
- allocate(psi_1h1p(N_int,2,n_det_1h1p))
+!  threshold_davidson = 1.d-12
- allocate(psi_1h2p(N_int,2,n_det_1h2p))
+!  soft_touch threshold_davidson davidson_criterion
- allocate(psi_ref_coef_out(N_det_generators,N_states))
+!  call diagonalize_CI
- allocate(psi_coef_1h1p(n_det_1h1p,N_states))
+!  call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p)
- allocate(psi_coef_1h2p(n_det_1h2p,N_states))
+!  allocate(psi_ref_out(N_int,2,N_det_generators))
- call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p)
+!  allocate(psi_1h1p(N_int,2,n_det_1h1p))
- call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
+!  allocate(psi_1h2p(N_int,2,n_det_1h2p))
-                                                        psi_1h1p,psi_coef_1h1p,n_det_1h1p)
+!  allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
- call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
+!  allocate(psi_ref_coef_out(N_det_generators,N_states))
-                                                        psi_1h2p,psi_coef_1h2p,n_det_1h2p)
+!  allocate(psi_coef_1h1p(n_det_1h1p,N_states))
-
+!  allocate(psi_coef_1h2p(n_det_1h2p,N_states))
- deallocate(psi_ref_out)
+!  allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
- deallocate(psi_1h1p)
+!  call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
- deallocate(psi_1h2p)
+!  call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
- deallocate(psi_ref_coef_out)
+!                                                         psi_1h1p,psi_coef_1h1p,n_det_1h1p)
- deallocate(psi_coef_1h1p)
+!  call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
- deallocate(psi_coef_1h2p)
+!                                                         psi_1h2p,psi_coef_1h2p,n_det_1h2p)
-
+!  call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
-end
+!                                                         psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
+! 
+!  deallocate(psi_ref_out)
+!  deallocate(psi_1h1p)
+!  deallocate(psi_1h2p)
+!  deallocate(psi_ref_coef_out)
+!  deallocate(psi_coef_1h1p)
+!  deallocate(psi_coef_1h2p)
+! 
+!  end
 
 

plugins/FOBOCI/collect_all_lmct.irp.f

@@ -0,0 +1,436 @@
+use bitmasks
+
+subroutine collect_lmct(hole_particle,n_couples)
+ implicit none
+ integer, intent(out) :: hole_particle(1000,2), n_couples
+ BEGIN_DOC
+ ! Collect all the couple holes/particles of the important LMCT
+ ! hole_particle(i,1) = ith hole
+ ! hole_particle(i,2) = ith particle
+ ! n_couples is the number of important excitations
+ END_DOC
+ print*,'COLLECTING THE PERTINENT LMCT (1h)'
+ double precision, allocatable :: tmp(:,:)
+ allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
+ tmp =  one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
+ integer :: i,j,iorb,jorb
+ n_couples = 0
+ do i = 1,n_act_orb
+  iorb = list_act(i)
+  do j = 1, n_inact_orb
+   jorb = list_inact(j)
+   if(dabs(tmp(iorb,jorb)).gt.1.d-2)then
+    n_couples +=1
+    hole_particle(n_couples,1) = jorb
+    hole_particle(n_couples,2) = iorb
+    print*,'DM'
+    print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
+   endif
+  enddo
+ enddo
+ deallocate(tmp)
+ print*,'number of meaning full couples of holes/particles '
+ print*,'n_couples = ',n_couples
+
+
+end
+
+
+subroutine collect_mlct(hole_particle,n_couples)
+ implicit none
+ integer, intent(out) :: hole_particle(1000,2), n_couples
+ BEGIN_DOC
+ ! Collect all the couple holes/particles of the important LMCT
+ ! hole_particle(i,1) = ith hole
+ ! hole_particle(i,2) = ith particle
+ ! n_couples is the number of important excitations
+ END_DOC
+ print*,'COLLECTING THE PERTINENT MLCT (1p)'
+ double precision, allocatable :: tmp(:,:)
+ allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
+ tmp =  one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
+ integer :: i,j,iorb,jorb
+ n_couples = 0
+ do i = 1,n_act_orb
+  iorb = list_act(i)
+  do j = 1, n_virt_orb
+   jorb = list_virt(j)
+   if(dabs(tmp(iorb,jorb)).gt.1.d-3)then
+    n_couples +=1
+    hole_particle(n_couples,1) = iorb
+    hole_particle(n_couples,2) = jorb
+    print*,'DM'
+    print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
+   endif
+  enddo
+ enddo
+ deallocate(tmp)
+ print*,'number of meaning full couples of holes/particles '
+ print*,'n_couples = ',n_couples
+
+
+end
+
+
+subroutine collect_lmct_mlct(hole_particle,n_couples)
+ implicit none
+ integer, intent(out) :: hole_particle(1000,2), n_couples
+ BEGIN_DOC
+ ! Collect all the couple holes/particles of the important LMCT
+ ! hole_particle(i,1) = ith hole
+ ! hole_particle(i,2) = ith particle
+ ! n_couples is the number of important excitations
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ print*,'COLLECTING THE PERTINENT LMCT (1h)'
+ print*,'AND THE PERTINENT MLCT (1p)'
+ allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
+ tmp =  one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
+ integer :: i,j,iorb,jorb
+ n_couples = 0
+ do i = 1,n_act_orb
+  iorb = list_act(i)
+  do j = 1, n_inact_orb
+   jorb = list_inact(j)
+   if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then
+    n_couples +=1
+    hole_particle(n_couples,1) = jorb
+    hole_particle(n_couples,2) = iorb
+    print*,'DM'
+    print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
+   endif
+  enddo
+  do j = 1, n_virt_orb
+   jorb = list_virt(j)
+   if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then
+    n_couples +=1
+    hole_particle(n_couples,1) = iorb
+    hole_particle(n_couples,2) = jorb
+    print*,'DM'
+    print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
+   endif
+  enddo
+ enddo
+ deallocate(tmp)
+ print*,'number of meaning full couples of holes/particles '
+ print*,'n_couples = ',n_couples
+
+
+end
+
+
+subroutine collect_1h1p(hole_particle,n_couples)
+ implicit none
+ integer, intent(out) :: hole_particle(1000,2), n_couples
+ BEGIN_DOC
+ ! Collect all the couple holes/particles of the important LMCT
+ ! hole_particle(i,1) = ith hole
+ ! hole_particle(i,2) = ith particle
+ ! n_couples is the number of important excitations
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ print*,'COLLECTING THE PERTINENT 1h1p'
+ allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
+ tmp =  one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
+ integer :: i,j,iorb,jorb
+ n_couples = 0
+ do i = 1,n_virt_orb
+  iorb = list_virt(i)
+  do j = 1, n_inact_orb
+   jorb = list_inact(j)
+   if(dabs(tmp(iorb,jorb)).gt.1.d-2)then
+    n_couples +=1
+    hole_particle(n_couples,1) = jorb
+    hole_particle(n_couples,2) = iorb
+    print*,'DM'
+    print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
+   endif
+  enddo
+ enddo
+ deallocate(tmp)
+ print*,'number of meaning full couples of holes/particles '
+ print*,'n_couples = ',n_couples
+
+
+end
+
+
+
+subroutine set_lmct_to_generators_restart 
+ implicit none
+ integer :: i,j,m,n,i_hole,i_particle
+ integer :: hole_particle(1000,2), n_couples
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: N_det_total,i_ok
+
+ call collect_lmct(hole_particle,n_couples)
+ call set_generators_to_generators_restart
+ N_det_total = N_det_generators_restart 
+ do i = 1, n_couples
+  i_hole     = hole_particle(i,1)
+  i_particle = hole_particle(i,2)
+  do m = 1, N_det_cas
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+   ! You excite the beta electron from i_hole to i_particle
+   print*,'i_hole,i_particle 2 = ',i_hole,i_particle 
+   call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
+   print*,'i_ok = ',i_ok
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+
+   ! You excite the alpha electron from i_hole to i_particle
+   print*,'i_hole,i_particle 1 = ',i_hole,i_particle 
+   call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
+   print*,'i_ok = ',i_ok
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+  enddo
+ enddo
+ N_det_generators = N_det_total
+ do i = 1, N_det_generators
+  psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
+ enddo
+ print*,'number of generators in total = ',N_det_generators
+ touch N_det_generators psi_coef_generators psi_det_generators
+end
+
+subroutine set_mlct_to_generators_restart 
+ implicit none
+ integer :: i,j,m,n,i_hole,i_particle
+ integer :: hole_particle(1000,2), n_couples
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: N_det_total,i_ok
+
+ call collect_mlct(hole_particle,n_couples)
+ call set_generators_to_generators_restart
+ N_det_total = N_det_generators_restart 
+ do i = 1, n_couples
+  i_hole     = hole_particle(i,1)
+  i_particle = hole_particle(i,2)
+  do m = 1, N_det_cas
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+   ! You excite the beta electron from i_hole to i_particle
+   print*,'i_hole,i_particle 2 = ',i_hole,i_particle 
+   call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
+   print*,'i_ok = ',i_ok
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+
+   ! You excite the alpha electron from i_hole to i_particle
+   print*,'i_hole,i_particle 1 = ',i_hole,i_particle 
+   call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
+   print*,'i_ok = ',i_ok
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+  enddo
+ enddo
+ N_det_generators = N_det_total
+ do i = 1, N_det_generators
+  psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
+ enddo
+ print*,'number of generators in total = ',N_det_generators
+ touch N_det_generators psi_coef_generators psi_det_generators
+end
+
+subroutine set_lmct_mlct_to_generators_restart 
+ implicit none
+ integer :: i,j,m,n,i_hole,i_particle
+ integer :: hole_particle(1000,2), n_couples
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: N_det_total,i_ok
+
+ call collect_lmct_mlct(hole_particle,n_couples)
+ call set_generators_to_generators_restart
+ N_det_total = N_det_generators_restart 
+ do i = 1, n_couples
+  i_hole     = hole_particle(i,1)
+  i_particle = hole_particle(i,2)
+  do m = 1, N_det_cas
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+   ! You excite the beta electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+
+   do n = 1, N_int
+    key_tmp(n,1) = psi_cas(n,1,m)
+    key_tmp(n,2) = psi_cas(n,2,m)
+   enddo
+
+   ! You excite the alpha electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
+     psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+  enddo
+ enddo
+ N_det_generators = N_det_total
+ do i = 1, N_det_generators
+  psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
+ enddo
+ print*,'number of generators in total = ',N_det_generators
+ touch N_det_generators psi_coef_generators psi_det_generators
+end
+
+subroutine set_lmct_mlct_to_psi_det
+ implicit none
+ integer :: i,j,m,n,i_hole,i_particle
+ integer :: hole_particle(1000,2), n_couples
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: N_det_total,i_ok
+
+ call collect_lmct_mlct(hole_particle,n_couples)
+ call set_psi_det_to_generators_restart
+ N_det_total = N_det_generators_restart 
+ do i = 1, n_couples
+  i_hole     = hole_particle(i,1)
+  i_particle = hole_particle(i,2)
+  do m = 1, N_det_generators_restart
+   do n = 1, N_int
+    key_tmp(n,1) = psi_det_generators_restart(n,1,m)
+    key_tmp(n,2) = psi_det_generators_restart(n,2,m)
+   enddo
+   ! You excite the beta electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det(n,1,N_det_total) = key_tmp(n,1)
+     psi_det(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+
+   do n = 1, N_int
+    key_tmp(n,1) = psi_det_generators_restart(n,1,m)
+    key_tmp(n,2) = psi_det_generators_restart(n,2,m)
+   enddo
+
+   ! You excite the alpha electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det(n,1,N_det_total) = key_tmp(n,1)
+     psi_det(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+  enddo
+ enddo
+ 
+ N_det = N_det_total
+ integer :: k
+ do k = 1, N_states
+  do i = 1, N_det
+   psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total))
+  enddo
+ enddo
+ SOFT_TOUCH N_det psi_det psi_coef
+ logical :: found_duplicates
+ call remove_duplicates_in_psi_det(found_duplicates)
+end
+
+subroutine set_1h1p_to_psi_det
+ implicit none
+ integer :: i,j,m,n,i_hole,i_particle
+ integer :: hole_particle(1000,2), n_couples
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: N_det_total,i_ok
+
+ call collect_1h1p(hole_particle,n_couples)
+ call set_psi_det_to_generators_restart
+ N_det_total = N_det_generators_restart 
+ do i = 1, n_couples
+  i_hole     = hole_particle(i,1)
+  i_particle = hole_particle(i,2)
+  do m = 1, N_det_generators_restart
+   do n = 1, N_int
+    key_tmp(n,1) = psi_det_generators_restart(n,1,m)
+    key_tmp(n,2) = psi_det_generators_restart(n,2,m)
+   enddo
+   ! You excite the beta electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det(n,1,N_det_total) = key_tmp(n,1)
+     psi_det(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+
+   do n = 1, N_int
+    key_tmp(n,1) = psi_det_generators_restart(n,1,m)
+    key_tmp(n,2) = psi_det_generators_restart(n,2,m)
+   enddo
+
+   ! You excite the alpha electron from i_hole to i_particle
+   call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
+   if(i_ok==1)then
+    N_det_total +=1
+    do n = 1, N_int
+     psi_det(n,1,N_det_total) = key_tmp(n,1)
+     psi_det(n,2,N_det_total) = key_tmp(n,2)
+    enddo
+   endif
+  enddo
+ enddo
+ 
+ N_det = N_det_total
+ integer :: k
+ do k = 1, N_states
+  do i = 1, N_det
+   psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total))
+  enddo
+ enddo
+ SOFT_TOUCH N_det psi_det psi_coef
+ logical :: found_duplicates
+ call remove_duplicates_in_psi_det(found_duplicates)
+end
+

plugins/FOBOCI/corr_energy_2h2p.irp.f

@@ -0,0 +1,425 @@
+ 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
+ integer :: i_hole,j_hole,k_part,l_part
+ double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
+ double precision :: diag_H_mat_elem
+ 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
+   k_part = list_virt(k)
+   do j = 1, n_inact_orb  ! alpha
+    j_hole = list_inact(j)
+    do l = 1, n_virt_orb  ! alpha
+     l_part = list_virt(l)
+
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     delta_e  = (ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
+
+     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
+ enddo
+
+ ! alpha alpha correlation energy
+ do i = 1, n_inact_orb
+  i_hole = list_inact(i)
+  do j = i+1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_virt_orb
+    k_part = list_virt(k)
+    do l = k+1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+     key_tmp = ref_bitmask
+     ispin = 1
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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))
+     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
+ enddo
+
+ ! beta beta correlation energy
+ do i = 1, n_inact_orb
+  i_hole = list_inact(i)
+  do j = i+1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_virt_orb
+    k_part = list_virt(k)
+    do l = k+1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 2
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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))
+     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
+ enddo
+
+END_PROVIDER 
+
+ 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]
+ use bitmasks
+ implicit none
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: i,j,k,l
+ integer :: i_hole,j_hole,k_part,l_part
+ double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
+ double precision :: diag_H_mat_elem
+ integer :: i_ok,ispin
+ ! Alpha - Beta correlation energy
+ total_corr_e_2h1p = 0.d0
+ 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
+   k_part = list_act(k)
+   do j = 1, n_inact_orb
+    j_hole = list_inact(j)
+    do l = 1, n_virt_orb
+     l_part = list_virt(l)
+
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     delta_e  = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
+
+     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 
+   enddo
+  enddo
+ enddo
+
+ ! Alpha Alpha spin correlation energy
+ do i = 1, n_inact_orb
+  i_hole = list_inact(i)
+  do j = i+1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_act_orb
+    k_part = list_act(k)
+    do l = 1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+     key_tmp = ref_bitmask
+     ispin = 1
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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))
+
+     total_corr_e_2h1p += contrib
+     corr_energy_2h1p_per_orb_aa(i_hole) += contrib
+     corr_energy_2h1p_per_orb_aa(l_part) += contrib
+    enddo 
+   enddo
+  enddo
+ enddo
+
+ ! Beta Beta correlation energy
+ do i = 1, n_inact_orb
+  i_hole = list_inact(i)
+  do j = i+1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_act_orb
+    k_part = list_act(k)
+    do l = 1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 2
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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
+     corr_energy_2h1p_per_orb_aa(l_part) += contrib
+    enddo 
+   enddo
+  enddo
+ enddo
+
+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]
+ use bitmasks
+ implicit none
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: i,j,k,l
+ integer :: i_hole,j_hole,k_part,l_part
+ double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
+ double precision :: diag_H_mat_elem
+ integer :: i_ok,ispin
+ ! Alpha - Beta correlation energy
+ total_corr_e_1h2p = 0.d0
+ corr_energy_1h2p_per_orb_ab = 0.d0
+ corr_energy_1h2p_per_orb_aa = 0.d0
+ corr_energy_1h2p_per_orb_bb = 0.d0
+ do i = 1, n_virt_orb
+  i_hole = list_virt(i)
+  do k = 1, n_act_orb
+   k_part = list_act(k)
+   do j = 1, n_inact_orb
+    j_hole = list_inact(j)
+    do l = 1, n_virt_orb
+     l_part = list_virt(l)
+
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     delta_e  = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
+
+     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_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
+ enddo
+
+ ! Alpha Alpha correlation energy
+ do i = 1, n_virt_orb
+  i_hole = list_virt(i)
+  do j = 1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_act_orb
+    k_part = list_act(k)
+    do l = i+1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+
+     key_tmp = ref_bitmask
+     ispin = 1
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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))
+     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
+ enddo
+
+ ! Beta Beta correlation energy
+ do i = 1, n_virt_orb
+  i_hole = list_virt(i)
+  do j = 1, n_inact_orb
+   j_hole = list_inact(j)
+   do k = 1, n_act_orb
+    k_part = list_act(k)
+    do l = i+1,n_virt_orb
+     l_part = list_virt(l)
+     hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
+     exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
+
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 2
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     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))
+     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
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, corr_energy_1h1p_spin_flip_per_orb, (mo_tot_num)]
+&BEGIN_PROVIDER [ double precision, total_corr_e_1h1p_spin_flip]
+ use bitmasks
+ implicit none
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: i,j,k,l
+ integer :: i_hole,j_hole,k_part,l_part
+ double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
+ double precision :: diag_H_mat_elem
+ integer :: i_ok,ispin
+ ! Alpha - Beta correlation energy
+ total_corr_e_1h1p_spin_flip = 0.d0
+ corr_energy_1h1p_spin_flip_per_orb = 0.d0
+ do i = 1, n_inact_orb
+  i_hole = list_inact(i)
+  do k = 1, n_act_orb
+   k_part = list_act(k)
+   do j = 1, n_act_orb
+    j_hole = list_act(j)
+    do l = 1, n_virt_orb
+     l_part = list_virt(l)
+
+     key_tmp = ref_bitmask
+     ispin = 2
+     call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     ispin = 1
+     call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
+     if(i_ok .ne.1)cycle
+     delta_e  = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
+
+     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_1h1p_spin_flip += contrib
+     corr_energy_1h1p_spin_flip_per_orb(i_hole) += contrib
+    enddo 
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 

plugins/FOBOCI/diag_fock_inactiv_virt.irp.f

@@ -3,6 +3,7 @@ subroutine diag_inactive_virt_and_update_mos
  integer :: i,j,i_inact,j_inact,i_virt,j_virt
  double precision :: tmp(mo_tot_num_align,mo_tot_num)
  character*(64) :: label
+ print*,'Diagonalizing the occ and virt Fock operator'
  tmp = 0.d0
  do i = 1, mo_tot_num
   tmp(i,i) = Fock_matrix_mo(i,i)
@@ -33,3 +34,50 @@ subroutine diag_inactive_virt_and_update_mos
 
 
 end
+
+subroutine diag_inactive_virt_new_and_update_mos
+ implicit none
+ integer :: i,j,i_inact,j_inact,i_virt,j_virt,k,k_act
+ double precision :: tmp(mo_tot_num_align,mo_tot_num),accu,get_mo_bielec_integral_schwartz
+ character*(64) :: label
+ tmp = 0.d0
+ do i = 1, mo_tot_num
+  tmp(i,i) = Fock_matrix_mo(i,i)
+ enddo
+ 
+ do i = 1, n_inact_orb
+  i_inact = list_inact(i)
+  do j = i+1, n_inact_orb 
+   j_inact = list_inact(j)
+   accu =0.d0
+   do k = 1, n_act_orb
+    k_act = list_act(k)
+    accu += get_mo_bielec_integral_schwartz(i_inact,k_act,j_inact,k_act,mo_integrals_map)
+    accu -= get_mo_bielec_integral_schwartz(i_inact,k_act,k_act,j_inact,mo_integrals_map)
+   enddo
+   tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact) + accu
+   tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact) + accu
+  enddo
+ enddo
+
+ do i = 1, n_virt_orb
+  i_virt = list_virt(i)
+  do j = i+1, n_virt_orb 
+   j_virt = list_virt(j)
+   accu =0.d0
+   do k = 1, n_act_orb
+    k_act = list_act(k)
+    accu += get_mo_bielec_integral_schwartz(i_virt,k_act,j_virt,k_act,mo_integrals_map)
+   enddo
+   tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt) - accu
+   tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt) - accu
+  enddo
+ enddo
+
+
+ label = "Canonical"
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
+ soft_touch mo_coef
+
+
+end

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
@@ -85,33 +85,6 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen
          delta_ij_generators_(idx(j), idx(k)) += contrib
        enddo 
      enddo
-!      H_matrix_tmp_bis(idx(k),idx(k)) += contrib
-!        H_matrix_tmp_bis(idx(k),idx(j)) += contrib
-!        H_matrix_tmp_bis(idx(j),idx(k)) += contrib
-!    do k = 1, Ndet_generators
-!     do j = 1, Ndet_generators
-!      H_matrix_tmp_bis(k,j) = H_matrix_tmp(k,j)
-!     enddo
-!    enddo
-!    double precision :: H_matrix_tmp_bis(Ndet_generators,Ndet_generators)
-!    double precision :: eigenvectors_bis(Ndet_generators,Ndet_generators), eigenvalues_bis(Ndet_generators)
-!    call lapack_diag(eigenvalues_bis,eigenvectors_bis,H_matrix_tmp_bis,Ndet_generators,Ndet_generators)
-!    print*,'f,lambda_i = ',f,lambda_i
-!    print*,'eigenvalues_bi(1)',eigenvalues_bis(1)
-!    print*,'eigenvalues      ',eigenvalues(1)
-!    do k = 1, Ndet_generators
-!     print*,'coef,coef_dres = ', eigenvectors(k,1), eigenvectors_bis(k,1)
-!    enddo
-!    pause
-!    accu = 0.d0
-!    do k = 1, Ndet_generators
-!     do j = 1, Ndet_generators
-!      accu += eigenvectors(k,1) * eigenvectors(j,1) * (H_matrix_tmp(k,j) + delta_ij_generators_(k,j))
-!     enddo
-!    enddo
-!    print*,'accu,eigv = ',accu,eigenvalues(1)
-!    pause
-
   enddo
 end
 

plugins/FOBOCI/fobo_scf.irp.f

@@ -0,0 +1,59 @@
+program foboscf
+ implicit none
+ call run_prepare
+ no_oa_or_av_opt = .True.
+ touch no_oa_or_av_opt
+ call routine_fobo_scf
+ call save_mos
+
+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
+
+subroutine routine_fobo_scf
+ implicit none
+ integer :: i,j
+ print*,''
+ print*,''
+ 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
+  call diag_inactive_virt_and_update_mos
+  call clear_mo_map
+  call provide_properties
+ enddo
+
+
+
+end

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,7 +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*,'--------------------------'
@@ -51,27 +55,91 @@ 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
-   ! so all the mono excitation on the new generators 
    allocate(dressing_matrix(N_det_generators,N_det_generators))
+   dressing_matrix = 0.d0
    if(.not.do_it_perturbative)then
-!   call all_single
+
-    dressing_matrix = 0.d0
     do k = 1, N_det_generators
      do l = 1, N_det_generators
        call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
        dressing_matrix(k,l) = hkl
      enddo
     enddo
-    double precision :: hkl
+    hkl = dressing_matrix(1,1)
-!   call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
+    do k = 1, N_det_generators
-!   call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
+      dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
-    call debug_det(reunion_of_bitmask,N_int)
+    enddo
+    print*,'Naked matrix'
+    do k = 1, N_det_generators
+     write(*,'(100(F12.5,X))')dressing_matrix(k,:)
+    enddo
+
+    ! Do all the single excitations on top of the CAS and 1h determinants
+    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 
+!   ! to introduce the Coulomb hole in the active space
+!   ! These are the 1h2p excitations that have the i_hole_osoci hole in common
+!   ! and the 2p if there is more than one electron in the active space
+!   do k = 1, N_int
+!    zero_bitmask(k,1) = 0_bit_kind
+!    zero_bitmask(k,2) = 0_bit_kind
+!   enddo
+!   ! hole is possible only in the orbital i_hole_osoci 
+!   call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
+!   call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
+!   ! and in the active space
+!   do k = 1, n_act_orb
+!    call set_bit_to_integer(list_act(k),zero_bitmask(1,1),N_int)
+!    call set_bit_to_integer(list_act(k),zero_bitmask(1,2),N_int)
+!   enddo
+!   call set_bitmask_hole_as_input(zero_bitmask)
+
+!   call set_bitmask_particl_as_input(reunion_of_bitmask)
+
+!   call all_1h2p
+!   call diagonalize_CI_SC2
+!   call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
+
+!   ! Change the mask of the holes and particles to perform all the 
+!   ! double excitations that from the orbital i_hole_osoci
+!   do k = 1, N_int
+!    zero_bitmask(k,1) = 0_bit_kind
+!    zero_bitmask(k,2) = 0_bit_kind
+!   enddo
+!   ! hole is possible only in the orbital i_hole_osoci 
+!   call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
+!   call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
+!   call set_bitmask_hole_as_input(zero_bitmask)
+
+!   call set_bitmask_particl_as_input(reunion_of_bitmask)
+
+!   call set_psi_det_to_generators
+!   call all_2h2p
+!   call diagonalize_CI_SC2
+    double precision :: hkl
+    call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
+    hkl = dressing_matrix(1,1)
+    do k = 1, N_det_generators
+      dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
+    enddo
+    print*,'Dressed matrix'
+    do k = 1, N_det_generators
+     write(*,'(100(F12.5,X))')dressing_matrix(k,:)
+    enddo
+!   call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
    endif
    call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci)
+
    do k = 1, N_states
     print*,'norm_tmp = ',norm_tmp(k)
     norm_total(k) += norm_tmp(k)
@@ -83,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
@@ -107,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))
@@ -122,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
@@ -132,24 +206,6 @@ subroutine FOBOCI_lmct_mlct_old_thr
    deallocate(dressing_matrix)
   enddo
  endif
-  if(.False.)then
-   print*,'LAST loop for all the 1h-1p'
-   print*,'--------------------------'
-   ! First set the current generators to the one of restart
-   call set_generators_to_generators_restart
-   call set_psi_det_to_generators
-   call initialize_bitmask_to_restart_ones
-   ! Impose that only the hole i_hole_osoci can be done
-   call set_bitmask_particl_as_input(inact_virt_bitmask)
-   call set_bitmask_hole_as_input(inact_virt_bitmask)
-!  call set_bitmask_particl_as_input(reunion_of_bitmask)
-!  call set_bitmask_hole_as_input(reunion_of_bitmask)
-   call all_single
-   call set_intermediate_normalization_1h1p(norm_tmp)
-   norm_total += norm_tmp
-   call update_density_matrix_osoci
-  endif
-
 
    print*,'norm_total = ',norm_total
    norm_total = norm_generators_restart
@@ -174,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))
@@ -216,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
@@ -250,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))
@@ -290,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/foboci_reunion.irp.f

@@ -0,0 +1,18 @@
+program osoci_program
+implicit none
+   do_it_perturbative = .True.
+   touch do_it_perturbative
+   call FOBOCI_lmct_mlct_old_thr
+   call provide_all_the_rest
+end
+subroutine provide_all_the_rest
+implicit none
+integer :: i
+   call update_one_body_dm_mo
+   call set_lmct_mlct_to_psi_det
+   call diagonalize_CI
+   call save_wavefunction
+
+
+end
+

plugins/FOBOCI/generators_restart_save.irp.f

@@ -1,126 +1,74 @@
-use bitmasks
 
+use bitmasks
+ 
 BEGIN_PROVIDER [ integer, N_det_generators_restart ]
  implicit none
  BEGIN_DOC
- ! Number of determinants in the wave function
+ ! Read the wave function 
  END_DOC
- logical                        :: exists
+ integer :: i
- character*64                   :: label
  integer, save :: ifirst = 0
-!if(ifirst == 0)then
+ double precision :: norm
- PROVIDE ezfio_filename
+ if(ifirst == 0)then
-   call ezfio_has_determinants_n_det(exists)
+  call ezfio_get_determinants_n_det(N_det_generators_restart)
-   print*,'exists = ',exists
-   if(.not.exists)then
-    print*,'The OSOCI needs a restart WF'
-    print*,'There are none in the EZFIO file ...'
-    print*,'Stopping ...'
-    stop
-   endif
-   print*,'passed N_det_generators_restart'
-   call ezfio_get_determinants_n_det(N_det_generators_restart)
- ASSERT (N_det_generators_restart > 0)
   ifirst = 1
-!endif
+ else
+  print*,'PB in generators_restart restart !!!'
+ endif
+ call write_int(output_determinants,N_det_generators_restart,'Number of generators_restart')
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,N_det_generators_restart) ]
 &BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ]
  implicit none
  BEGIN_DOC
- ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
+ ! read wf
- ! is empty
+ ! 
  END_DOC
-  integer                        :: i
+ integer                        :: i, k
-  logical                        :: exists
-  character*64                   :: label
-  
  integer, save :: ifirst = 0
-!if(ifirst == 0)then
+ double precision, allocatable  :: psi_coef_read(:,:)
-  provide N_det_generators_restart
+ if(ifirst == 0)then
-  if(.True.)then
+  call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
-    call ezfio_has_determinants_N_int(exists)
+   do k = 1, N_int
-    if (exists) then
+    ref_generators_restart(k,1) = psi_det_generators_restart(k,1,1)
-     call ezfio_has_determinants_bit_kind(exists)
+    ref_generators_restart(k,2) = psi_det_generators_restart(k,2,1)
-     if (exists) then
+   enddo
-      call ezfio_has_determinants_N_det(exists)
+   allocate (psi_coef_read(N_det_generators_restart,N_states))
-      if (exists) then
+   call ezfio_get_determinants_psi_coef(psi_coef_read)
-       call ezfio_has_determinants_N_states(exists)
+   do k = 1, N_states
-       if (exists) then
+    do i = 1, N_det_generators_restart
-        call ezfio_has_determinants_psi_det(exists)
+     psi_coef_generators_restart(i,k) = psi_coef_read(i,k)
-       endif
+    enddo
-      endif
+   enddo
-     endif
-    endif
-
-    if(.not.exists)then
-     print*,'The OSOCI needs a restart WF'
-     print*,'There are none in the EZFIO file ...'
-     print*,'Stopping ...'
-     stop
-    endif
-    print*,'passed psi_det_generators_restart'
-   
-      call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
-      do i = 1, N_int
-        ref_generators_restart(i,1) = psi_det_generators_restart(i,1,1)
-        ref_generators_restart(i,2) = psi_det_generators_restart(i,2,1)
-      enddo
-  endif
   ifirst = 1
-!endif
+  deallocate(psi_coef_read)
+ else 
+  print*,'PB in generators_restart restart !!!'
+ endif
 
 END_PROVIDER
 
-
+BEGIN_PROVIDER [ integer, size_select_max]
-
+ implicit none
-BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (psi_det_size,N_states_diag) ]
+ BEGIN_DOC
-  implicit none
+ ! Size of the select_max array
-  BEGIN_DOC
+ END_DOC
-  ! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
+ size_select_max = 10000
-  ! is empty
-  END_DOC
-  
-  integer                        :: i,k, N_int2
-  logical                        :: exists
-  double precision, allocatable  :: psi_coef_read(:,:)
-  character*(64)                 :: label
-
- integer, save :: ifirst = 0
-!if(ifirst == 0)then
-  psi_coef_generators_restart = 0.d0
-  do i=1,N_states_diag
-    psi_coef_generators_restart(i,i) = 1.d0
-  enddo
-
-    call ezfio_has_determinants_psi_coef(exists)
-    
-   if(.not.exists)then
-    print*,'The OSOCI needs a restart WF'
-    print*,'There are none in the EZFIO file ...'
-    print*,'Stopping ...'
-    stop
-   endif
-   print*,'passed  psi_coef_generators_restart'
-    
-    if (exists) then
-      
-      allocate (psi_coef_read(N_det_generators_restart,N_states))
-      call ezfio_get_determinants_psi_coef(psi_coef_read)
-      do k=1,N_states
-        do i=1,N_det_generators_restart
-          psi_coef_generators_restart(i,k) = psi_coef_read(i,k)
-        enddo
-      enddo
-      deallocate(psi_coef_read)
-      
-    endif
-  ifirst = 1
-!endif
-    
-    
-  
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
+ implicit none
+ BEGIN_DOC
+ ! Memo to skip useless selectors
+ END_DOC
+ select_max = huge(1.d0)
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, N_det_generators ]
+&BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,10000) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (10000,N_states) ]
+
+END_PROVIDER

plugins/FOBOCI/hcc_1h1p.irp.f

@@ -0,0 +1,83 @@
+program test_sc2
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+
+
+end
+
+subroutine routine
+  implicit none
+  double precision, allocatable :: energies(:),diag_H_elements(:)
+  double precision, allocatable :: H_matrix(:,:)
+  allocate(energies(N_states),diag_H_elements(N_det))
+  call diagonalize_CI
+  call test_hcc
+  call test_mulliken
+! call SC2_1h1p(psi_det,psi_coef,energies, &
+!       diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+  allocate(H_matrix(N_det,N_det))
+  call SC2_1h1p_full(psi_det,psi_coef,energies, &
+        H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+  deallocate(H_matrix)
+  integer :: i,j
+  double precision :: accu,coef_hf
+! coef_hf = 1.d0/psi_coef(1,1)
+! do i = 1, N_det
+!  psi_coef(i,1) *= coef_hf
+! enddo
+  touch psi_coef
+  call pouet
+end
+
+subroutine pouet
+  implicit none
+  double precision :: accu,coef_hf
+! provide one_body_dm_mo_alpha one_body_dm_mo_beta
+! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int)
+! touch one_body_dm_mo_alpha one_body_dm_mo_beta 
+  call test_hcc
+  call test_mulliken
+! call save_wavefunction
+
+end
+
+subroutine test_hcc
+ implicit none
+ double precision :: accu
+ integer :: i,j
+ print*,'Z               AU           GAUSS              MHZ             cm^-1'
+ do i = 1, nucl_num
+  write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
+ enddo
+
+end
+
+subroutine test_mulliken
+ double precision :: accu
+ integer :: i
+ integer :: j
+ accu= 0.d0
+ do i = 1, nucl_num
+  print*,i,nucl_charge(i),mulliken_spin_densities(i)
+  accu += mulliken_spin_densities(i)
+ enddo
+ print*,'Sum of Mulliken SD = ',accu
+!print*,'AO SPIN POPULATIONS'
+ accu = 0.d0
+!do i = 1, ao_num
+! accu += spin_gross_orbital_product(i)
+! write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
+!enddo
+!print*,'sum = ',accu
+!accu = 0.d0
+!print*,'Angular momentum analysis'
+!do i = 0,  ao_l_max
+! accu += spin_population_angular_momentum(i)
+! print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
+!print*,'sum = ',accu
+!enddo
+
+end
+

plugins/FOBOCI/modify_generators.irp.f

@@ -6,6 +6,7 @@ subroutine set_generators_to_psi_det
  END_DOC
  N_det_generators = N_det
  integer :: i,k 
+ print*,'N_det = ',N_det
  do i=1,N_det_generators
    do k=1,N_int
      psi_det_generators(k,1,i) = psi_det(k,1,i)

plugins/FOBOCI/new_approach.irp.f

@@ -24,6 +24,7 @@ subroutine new_approach
  double precision, allocatable :: dressing_matrix_1h1p(:,:)
  double precision, allocatable :: dressing_matrix_2h1p(:,:)
  double precision, allocatable :: dressing_matrix_1h2p(:,:)
+ double precision, allocatable :: dressing_matrix_extra_1h_or_1p(:,:)
  double precision, allocatable :: H_matrix_tmp(:,:)
  logical :: verbose,is_ok
 
@@ -45,7 +46,7 @@ subroutine new_approach
 
 
  verbose = .True.
- threshold = threshold_singles
+ threshold = threshold_lmct
  print*,'threshold = ',threshold
  thr = 1.d-12
  print*,''
@@ -81,12 +82,14 @@ subroutine new_approach
    ! so all the mono excitation on the new generators 
    allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators))
    allocate(dressing_matrix_2h1p(N_det_generators,N_det_generators))
+   allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators))
    dressing_matrix_1h1p = 0.d0
    dressing_matrix_2h1p = 0.d0
+   dressing_matrix_extra_1h_or_1p = 0.d0
    if(.not.do_it_perturbative)then
     n_good_hole +=1
 !   call all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
-    call all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
+    call all_single_for_1h(i_hole_foboci,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p)
     allocate(H_matrix_tmp(N_det_generators,N_det_generators))
     do j = 1,N_det_generators
      do k = 1, N_det_generators
@@ -96,7 +99,7 @@ subroutine new_approach
     enddo
     do j = 1, N_det_generators
      do k = 1, N_det_generators
-      H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k)
+      H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + dressing_matrix_extra_1h_or_1p(j,k)
      enddo
     enddo
     hjk = H_matrix_tmp(1,1)
@@ -130,6 +133,7 @@ subroutine new_approach
    endif
    deallocate(dressing_matrix_1h1p)
    deallocate(dressing_matrix_2h1p)
+   deallocate(dressing_matrix_extra_1h_or_1p)
  enddo
  
  print*,''
@@ -155,12 +159,14 @@ subroutine new_approach
    ! so all the mono excitation on the new generators 
    allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators))
    allocate(dressing_matrix_1h2p(N_det_generators,N_det_generators))
+   allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators))
    dressing_matrix_1h1p = 0.d0
    dressing_matrix_1h2p = 0.d0
+   dressing_matrix_extra_1h_or_1p = 0.d0
    if(.not.do_it_perturbative)then
     n_good_hole +=1
 !   call all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
-    call all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
+    call all_single_for_1p(i_particl_osoci,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
     allocate(H_matrix_tmp(N_det_generators,N_det_generators))
     do j = 1,N_det_generators
      do k = 1, N_det_generators
@@ -170,7 +176,7 @@ subroutine new_approach
     enddo
     do j = 1, N_det_generators
      do k = 1, N_det_generators
-      H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k)
+      H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + dressing_matrix_extra_1h_or_1p(j,k)
      enddo
     enddo
     hjk = H_matrix_tmp(1,1)
@@ -205,7 +211,10 @@ subroutine new_approach
    endif
    deallocate(dressing_matrix_1h1p)
    deallocate(dressing_matrix_1h2p)
+   deallocate(dressing_matrix_extra_1h_or_1p)
  enddo
+
+
  double precision, allocatable :: H_matrix_total(:,:)
  integer :: n_det_total
  n_det_total = N_det_generators_restart + n_good_det
@@ -221,7 +230,7 @@ subroutine new_approach
    !!! Adding the averaged dressing coming from the 1h1p that are redundant for each of the "n_good_hole" 1h
    H_matrix_total(i,j) += dressing_matrix_restart_1h1p(i,j)/dble(n_good_hole+n_good_particl)  
    !!! Adding the dressing coming from the 2h1p that are not redundant for the any of CI calculations
-   H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j)
+   H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + dressing_matrix_restart_1h2p(i,j)
   enddo
  enddo
  do i = 1, n_good_det
@@ -244,25 +253,79 @@ subroutine new_approach
    H_matrix_total(n_det_generators_restart+j,n_det_generators_restart+i) = hij 
   enddo
  enddo
- print*,'H matrix to diagonalize'
+
- double precision :: href
+ ! Adding the correlation energy 
- href = H_matrix_total(1,1)
+ logical :: orb_taken_good_det(mo_tot_num)
- do i = 1, n_det_total
+ double precision :: phase
-  H_matrix_total(i,i) -= href
+ integer :: n_h,n_p,number_of_holes,number_of_particles
+ integer :: exc(0:2,2,2)
+ integer :: degree
+ integer :: h1,h2,p1,p2,s1,s2
+ logical, allocatable :: one_hole_or_one_p(:)
+ integer, allocatable :: holes_or_particle(:)
+ allocate(one_hole_or_one_p(n_good_det), holes_or_particle(n_good_det))
+ orb_taken_good_det = .False.
+ do i = 1, n_good_det
+  n_h = number_of_holes(psi_good_det(1,1,i))
+  n_p = number_of_particles(psi_good_det(1,1,i))
+  call get_excitation(ref_bitmask,psi_good_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 == 1)then
+   orb_taken_good_det(h1) = .True.
+   one_hole_or_one_p(i) = .True.
+   holes_or_particle(i) = h1
+  endif
+  if(n_h == 1 .and. n_p == 0)then
+   orb_taken_good_det(p1) = .True.
+   one_hole_or_one_p(i) = .False.
+   holes_or_particle(i) = p1
+  endif
  enddo
- do i = 1, n_det_total
+ 
-   write(*,'(100(X,F16.8))')H_matrix_total(i,:)
+  do i = 1, N_det_generators_restart
- enddo
+   ! Add the 2h2p, 2h1p and 1h2p correlation energy
- double precision, allocatable :: eigvalues(:),eigvectors(:,:)
+   H_matrix_total(i,i) += total_corr_e_2h2p + total_corr_e_2h1p + total_corr_e_1h2p + total_corr_e_1h1p_spin_flip
- allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total))
+   ! Substract the 2h1p part that have already been taken into account
- call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total)
+   do j = 1, n_inact_orb
- print*,'e_dressed  = ',eigvalues(1) + nuclear_repulsion + href
+    iorb = list_inact(j)
- do i = 1, n_det_total
+    if(.not.orb_taken_good_det(iorb))cycle
-  print*,'coef = ',eigvectors(i,1)
+    H_matrix_total(i,i) -= corr_energy_2h1p_per_orb_ab(iorb) - corr_energy_2h1p_per_orb_bb(iorb) - corr_energy_1h1p_spin_flip_per_orb(iorb)
- enddo
+   enddo
- integer(bit_kind), allocatable :: psi_det_final(:,:,:)
+   ! Substract the 1h2p part that have already been taken into account
- double precision, allocatable :: psi_coef_final(:,:)
+   do j = 1, n_virt_orb
- double precision :: norm
+    iorb = list_virt(j)
+    if(.not.orb_taken_good_det(iorb))cycle
+    H_matrix_total(i,i) -= corr_energy_1h2p_per_orb_ab(iorb) - corr_energy_1h2p_per_orb_aa(iorb)
+   enddo
+  enddo
+  
+  do i = 1, N_good_det
+   ! Repeat the 2h2p correlation energy
+   H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_2h2p 
+   ! Substract the part that can not be repeated
+   ! If it is a 1h
+   if(one_hole_or_one_p(i))then
+    ! 2h2p 
+    H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) &  
+                                                                             -corr_energy_2h2p_per_orb_bb(holes_or_particle(i))  
+    ! You can repeat a certain part of the 1h2p correlation energy
+    ! that is everything except the part that involves the hole of the 1h
+    H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_1h2p 
+    H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_1h2p_per_orb_ab(holes_or_particle(i)) &  
+                                                                             -corr_energy_1h2p_per_orb_bb(holes_or_particle(i))  
+   
+   else
+    ! 2h2p 
+    H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) &  
+                                                                             -corr_energy_2h2p_per_orb_aa(holes_or_particle(i))  
+    ! You can repeat a certain part of the 2h1p correlation energy
+    ! that is everything except the part that involves the hole of the 1p
+    ! 2h1p
+    H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h1p_per_orb_ab(holes_or_particle(i)) &  
+                                                                             -corr_energy_2h1p_per_orb_aa(holes_or_particle(i))  
+   endif
+  enddo
+
  allocate(psi_coef_final(n_det_total, N_states))
  allocate(psi_det_final(N_int,2,n_det_total))
  do i = 1, N_det_generators_restart
@@ -277,22 +340,222 @@ subroutine new_approach
    psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i)
   enddo
  enddo
- norm = 0.d0
+
+
+ double precision :: href
+ double precision, allocatable :: eigvalues(:),eigvectors(:,:)
+ integer(bit_kind), allocatable :: psi_det_final(:,:,:)
+ double precision, allocatable :: psi_coef_final(:,:)
+ double precision :: norm
+ allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total))
+
+ call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total)
+ print*,''
+ print*,''
+ print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1)
+ print*,'e_dressed  = ',eigvalues(1) + nuclear_repulsion
  do i = 1, n_det_total
-  do j = 1, N_states
+  print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1)
-   psi_coef_final(i,j) = eigvectors(i,j)
-  enddo
-  norm += psi_coef_final(i,1)**2
-! call debug_det(psi_det_final(1, 1, i), N_int)
  enddo
- print*,'norm = ',norm
+
+ integer(bit_kind), allocatable :: psi_det_remaining_1h_or_1p(:,:,:)
+ integer(bit_kind), allocatable :: key_tmp(:,:)
+ integer :: n_det_remaining_1h_or_1p
+ integer :: ispin,i_ok
+ allocate(key_tmp(N_int,2),psi_det_remaining_1h_or_1p(N_int,2,n_inact_orb*n_act_orb+n_virt_orb*n_act_orb))
+ logical :: is_already_present
+ logical, allocatable :: one_hole_or_one_p_bis(:)
+ integer, allocatable :: holes_or_particle_bis(:)
+ double precision,allocatable :: H_array(:)
+ allocate(one_hole_or_one_p_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb), holes_or_particle_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb))
+ allocate(H_array(n_det_total))
+ ! Dressing with the remaining 1h determinants
+ print*,''
+ print*,''
+ print*,'Dressing with the remaining 1h determinants'
+ n_det_remaining_1h_or_1p = 0
+ do i = 1, n_inact_orb
+  iorb = list_inact(i)
+  if(orb_taken_good_det(iorb))cycle
+  do j = 1, n_act_orb
+   jorb = list_act(j)
+   ispin = 2
+   key_tmp = ref_bitmask
+   call do_mono_excitation(key_tmp,iorb,jorb,ispin,i_ok)
+   if(i_ok .ne.1)cycle
+   is_already_present = .False.
+   H_array = 0.d0
+   call i_h_j(key_tmp,key_tmp,N_int,hij)
+   href = ref_bitmask_energy - hij
+   href = 1.d0/href
+   do k = 1, n_det_total
+    call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int)
+    if(degree == 0)then
+     is_already_present = .True.
+     exit
+    endif
+   enddo
+   if(is_already_present)cycle
+   n_det_remaining_1h_or_1p +=1
+   one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .True.
+   holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb
+   do k = 1, N_int
+    psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1)
+    psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2)
+   enddo
+ ! do k = 1, n_det_total
+ !   call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij)
+ !   H_array(k) = hij
+ ! enddo
+ ! do k = 1, n_det_total 
+ !  do l = 1, n_det_total
+ !   H_matrix_total(k,l) += H_array(k) * H_array(l) * href
+ !  enddo
+ ! enddo
+  enddo
+ enddo
+ ! Dressing with the remaining 1p determinants
+ print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p
+ print*,'Dressing with the remaining 1p determinants'
+ do i = 1, n_virt_orb
+  iorb = list_virt(i)
+  if(orb_taken_good_det(iorb))cycle
+  do j = 1, n_act_orb
+   jorb = list_act(j)
+   ispin = 1
+   key_tmp = ref_bitmask
+   call do_mono_excitation(key_tmp,jorb,iorb,ispin,i_ok)
+   if(i_ok .ne.1)cycle
+   is_already_present = .False.
+   H_array = 0.d0
+   call i_h_j(key_tmp,key_tmp,N_int,hij)
+   href = ref_bitmask_energy - hij
+   href = 1.d0/href
+   do k = 1, n_det_total
+    call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int)
+    if(degree == 0)then
+     is_already_present = .True.
+     exit
+    endif
+   enddo
+   if(is_already_present)cycle
+   n_det_remaining_1h_or_1p +=1
+   one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .False.
+   holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb
+   do k = 1, N_int
+    psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1)
+    psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2)
+   enddo
+!  do k = 1, n_det_total
+!    call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij)
+!    H_array(k) = hij
+!  enddo
+!  do k = 1, n_det_total 
+!   do l = 1, n_det_total
+!    H_matrix_total(k,l) += H_array(k) * H_array(l) * href
+!   enddo
+!  enddo
+  enddo
+ enddo
+ print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p
+ deallocate(key_tmp,H_array)
+
+ double precision, allocatable :: eigvalues_bis(:),eigvectors_bis(:,:),H_matrix_total_bis(:,:)
+ integer :: n_det_final
+ n_det_final = n_det_total + n_det_remaining_1h_or_1p
+ allocate(eigvalues_bis(n_det_final),eigvectors_bis(n_det_final,n_det_final),H_matrix_total_bis(n_det_final,n_det_final))
+ print*,'passed the allocate, building the big matrix'
+ do i = 1, n_det_total 
+  do j = 1, n_det_total
+   H_matrix_total_bis(i,j) = H_matrix_total(i,j)
+  enddo
+ enddo
+ do i = 1, n_det_remaining_1h_or_1p
+  do j = 1, n_det_remaining_1h_or_1p
+   call i_h_j(psi_det_remaining_1h_or_1p(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij)
+   H_matrix_total_bis(n_det_total+i,n_det_total+j) = hij
+  enddo
+ enddo
+ do i = 1, n_det_total
+  do j = 1, n_det_remaining_1h_or_1p
+   call i_h_j(psi_det_final(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij)
+   H_matrix_total_bis(i,n_det_total+j) = hij 
+   H_matrix_total_bis(n_det_total+j,i) = hij 
+  enddo
+ enddo
+ print*,'passed the matrix'
+ do i = 1, n_det_remaining_1h_or_1p
+   if(one_hole_or_one_p_bis(i))then
+    H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) &
+                                                                         -corr_energy_2h2p_per_orb_bb(holes_or_particle_bis(i))  
+    H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_1h2p_per_orb_ab(holes_or_particle_bis(i)) &  
+                                                                         -corr_energy_1h2p_per_orb_bb(holes_or_particle_bis(i))  
+   else
+    H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) &
+                                                                         -corr_energy_2h2p_per_orb_aa(holes_or_particle_bis(i))  
+    H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_2h1p_per_orb_ab(holes_or_particle_bis(i)) &  
+                                                                         -corr_energy_2h1p_per_orb_aa(holes_or_particle_bis(i))  
+ 
+   endif
+ enddo
+ do i = 2, n_det_final
+  do  j = i+1, n_det_final
+   H_matrix_total_bis(i,j) = 0.d0
+   H_matrix_total_bis(j,i) = 0.d0
+  enddo
+ enddo
+ do i = 1, n_det_final
+   write(*,'(500(F10.5,X))')H_matrix_total_bis(i,:)
+ enddo
+ call lapack_diag(eigvalues_bis,eigvectors_bis,H_matrix_total_bis,n_det_final,n_det_final)
+ print*,'e_dressed  = ',eigvalues_bis(1) + nuclear_repulsion
+ do i = 1, n_det_final
+  print*,'coef = ',eigvectors_bis(i,1),H_matrix_total_bis(i,i) - H_matrix_total_bis(1,1)
+ enddo
+ do j = 1, N_states
+  do i = 1, n_det_total
+   psi_coef_final(i,j) = eigvectors_bis(i,j)
+   norm += psi_coef_final(i,j)**2
+  enddo
+  norm = 1.d0/dsqrt(norm)
+  do i = 1, n_det_total
+   psi_coef_final(i,j) = psi_coef_final(i,j) * norm
+  enddo
+ enddo
+
+
+ deallocate(eigvalues_bis,eigvectors_bis,H_matrix_total_bis)
+
+
+!print*,'H matrix to diagonalize'
+!href = H_matrix_total(1,1)
+!do i = 1, n_det_total
+! H_matrix_total(i,i) -= href
+!enddo
+!do i = 1, n_det_total
+!  write(*,'(100(X,F16.8))')H_matrix_total(i,:)
+!enddo
+!call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total)
+!print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1)
+!print*,'e_dressed  = ',eigvalues(1) + nuclear_repulsion
+!do i = 1, n_det_total
+! print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1)
+!enddo
+!norm = 0.d0
+!do i = 1, n_det_total
+! do j = 1, N_states
+!  psi_coef_final(i,j) = eigvectors(i,j)
+! enddo
+! norm += psi_coef_final(i,1)**2
+!enddo
+!print*,'norm = ',norm
   
  call set_psi_det_as_input_psi(n_det_total,psi_det_final,psi_coef_final)
- print*,''
+
-!do i = 1, N_det
+ do i = 1, N_det
-!  call debug_det(psi_det(1,1,i),N_int)
+   call debug_det(psi_det(1,1,i),N_int)
-!  print*,'coef = ',psi_coef(i,1)
+   print*,'coef = ',psi_coef(i,1)
-!enddo
+ enddo
  provide one_body_dm_mo
 
  integer :: i_core,iorb,jorb,i_inact,j_inact,i_virt,j_virt,j_core
@@ -360,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/new_new_approach.irp.f

@@ -0,0 +1,132 @@
+program test_new_new
+ implicit none 
+ read_wf = .True.
+ touch read_wf
+ call test
+end
+
+
+subroutine test
+ implicit none
+ integer :: i,j,k,l
+ call diagonalize_CI 
+ call set_generators_to_psi_det
+ print*,'Initial coefficients'
+ do i = 1, N_det
+  print*,''
+  call debug_det(psi_det(1,1,i),N_int)
+  print*,'psi_coef = ',psi_coef(i,1) 
+  print*,''
+ enddo
+ double precision, allocatable :: dressing_matrix(:,:)
+ double precision :: hij 
+ double precision :: phase
+ integer :: n_h,n_p,number_of_holes,number_of_particles
+ integer :: exc(0:2,2,2)
+ integer :: degree
+ integer :: h1,h2,p1,p2,s1,s2
+ allocate(dressing_matrix(N_det_generators,N_det_generators))
+ do i = 1, N_det_generators
+  do j = 1, N_det_generators 
+   call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij)
+   dressing_matrix(i,j) = hij
+  enddo
+ enddo
+ href =  dressing_matrix(1,1)
+ print*,'Diagonal part of the dressing'
+ do i = 1, N_det_generators
+  print*,'delta e = ',dressing_matrix(i,i) - href
+ enddo
+ call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
+ double precision :: href
+ print*,''
+ ! One considers that the following excitation classes are not repeatable on the 1h and 1p determinants : 
+ !   + 1h1p spin flip
+ !   + 2h1p 
+ !   + 1h2p 
+ ! But the 2h2p are correctly taken into account 
+!dressing_matrix(1,1) += total_corr_e_1h2p + total_corr_e_2h1p + total_corr_e_1h1p_spin_flip
+!do i = 1, N_det_generators
+! dressing_matrix(i,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))
+!  if(n_h == 1 .and. n_p ==0)then
+!   
+!   call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int)
+!   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+!   print*,''
+!   print*,' 1h det '
+!   print*,''
+!   call debug_det(psi_det_generators(1,1,i),N_int)
+!   print*,'h1,p1 = ',h1,p1
+!   print*,'total_corr_e_2h2p         ',total_corr_e_2h2p
+!   print*,'corr_energy_2h2p_per_orb_ab(h1)',corr_energy_2h2p_per_orb_ab(h1) 
+!   print*,'corr_energy_2h2p_per_orb_bb(h1)',corr_energy_2h2p_per_orb_bb(h1) 
+!   dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1)
+!   dressing_matrix(1,1) += -corr_energy_2h1p_per_orb_aa(h1) - corr_energy_2h1p_per_orb_ab(h1) -corr_energy_2h1p_per_orb_bb(h1) & 
+!                           -corr_energy_1h1p_spin_flip_per_orb(h1)
+!  endif
+!  if(n_h == 0 .and. n_p ==1)then
+!   call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int)
+!   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+!   print*,''
+!   print*,' 1p det '
+!   print*,''
+!   call debug_det(psi_det_generators(1,1,i),N_int)
+!   print*,'h1,p1 = ',h1,p1
+!   print*,'total_corr_e_2h2p         ',total_corr_e_2h2p
+!   print*,'corr_energy_2h2p_per_orb_ab(p1)',corr_energy_2h2p_per_orb_ab(p1) 
+!   print*,'corr_energy_2h2p_per_orb_aa(p1)',corr_energy_2h2p_per_orb_aa(p1) 
+!   dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(p1) - corr_energy_2h2p_per_orb_aa(p1)
+!   dressing_matrix(1,1) += -corr_energy_1h2p_per_orb_aa(p1) - corr_energy_1h2p_per_orb_ab(p1) -corr_energy_1h2p_per_orb_bb(p1)
+!  endif
+!enddo
+!href =  dressing_matrix(1,1)
+!print*,'Diagonal part of the dressing'
+!do i = 1, N_det_generators
+! print*,'delta e = ',dressing_matrix(i,i) - href
+!enddo
+ call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
+ print*,'After dressing matrix'
+ print*,''
+ print*,''
+ do i = 1, N_det
+  print*,'psi_coef = ',psi_coef(i,1) 
+ enddo
+!print*,''
+!print*,''
+!print*,'Canceling the dressing part of the interaction between 1h and 1p'
+!do i = 2, N_det_generators
+! do j = i+1, N_det_generators
+!  call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij)
+!  dressing_matrix(i,j) = hij
+!  dressing_matrix(j,i) = hij 
+! enddo
+!enddo
+!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
+!print*,''
+!print*,''
+!do i = 1, N_det
+! print*,'psi_coef = ',psi_coef(i,1) 
+!enddo
+!print*,''
+!print*,''
+!print*,'Canceling the interaction between 1h and 1p'
+
+!print*,''
+!print*,''
+!do i = 2, N_det_generators
+! do j = i+1, N_det_generators
+!  dressing_matrix(i,j) = 0.d0
+!  dressing_matrix(j,i) = 0.d0
+! enddo
+!enddo
+!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
+!do i = 1, N_det
+! print*,'psi_coef = ',psi_coef(i,1) 
+!enddo
+ call save_natural_mos
+ deallocate(dressing_matrix)
+
+
+end

plugins/FOBOCI/routines_dressing.irp.f

@@ -55,15 +55,11 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det
    i_pert = 0
   endif
   do j = 1, ndet_generators_input
-   if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
+   if(dabs(H_array(j)*lambda_i).gt.0.1d0)then
     i_pert = 1
     exit
    endif
   enddo
-! print*,''
-! print*,'lambda_i,f = ',lambda_i,f
-! print*,'i_pert = ',i_pert
-! print*,''
   if(i_pert==1)then
    lambda_i = f
    i_pert_count +=1
@@ -79,9 +75,122 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det
     enddo 
   enddo
  enddo
+ href = dressing_matrix(1,1)
+ print*,'Diagonal part of the dressing'
+ do i = 1, ndet_generators_input
+  print*,'delta e = ',dressing_matrix(i,i) - href
+ enddo
 !print*,'i_pert_count = ',i_pert_count
 end
 
+subroutine update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,H_jj_in)
+ use bitmasks
+ implicit none
+ integer, intent(in) :: ndet_generators_input
+ integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
+ double precision, intent(in)  :: H_jj_in(N_det)
+ double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
+ integer :: i,j,n_det_ref_tmp,degree
+ double precision :: href
+ n_det_ref_tmp = 0
+ do i = 1, N_det
+  do j = 1, Ndet_generators_input
+   call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int)  
+   if(degree == 0)then
+    dressing_matrix(j,j) += H_jj_in(i)
+    n_det_ref_tmp +=1
+    exit
+   endif
+  enddo
+ enddo
+ if( ndet_generators_input .ne. n_det_ref_tmp)then
+  print*,'Problem !!!! '
+  print*,' ndet_generators .ne. n_det_ref_tmp !!!'
+  print*,'ndet_generators,n_det_ref_tmp'
+  print*,ndet_generators_input,n_det_ref_tmp 
+  stop
+ endif
+ 
+ href = dressing_matrix(1,1)
+ print*,''
+ print*,'Update with the SC2 dressing'
+ print*,''
+ print*,'Diagonal part of the dressing'
+ do i = 1, ndet_generators_input
+  print*,'delta e = ',dressing_matrix(i,i) - href
+ enddo
+end
+
+subroutine provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, &
+                                                           psi_det_outer_input,psi_coef_outer_input,n_det_outer_input)
+ use bitmasks
+ implicit none
+ integer, intent(in) :: n_det_ref_input
+ integer(bit_kind), intent(in) :: psi_det_ref_input(N_int,2,n_det_ref_input)
+ double precision, intent(in) :: psi_coef_ref_input(n_det_ref_input,N_states)
+ integer, intent(in) :: n_det_outer_input
+ integer(bit_kind), intent(in) :: psi_det_outer_input(N_int,2,n_det_outer_input)
+ double precision, intent(in) :: psi_coef_outer_input(n_det_outer_input,N_states)
+
+ double precision, intent(inout) :: dressing_matrix(n_det_ref_input,n_det_ref_input)
+
+
+ integer :: i_pert, i_pert_count,i,j,k
+ double precision :: f,href,hka,lambda_i
+ double precision :: H_array(n_det_ref_input),accu
+ integer :: n_h_out,n_p_out,n_p_in,n_h_in,number_of_holes,number_of_particles
+ call i_h_j(psi_det_ref_input(1,1,1),psi_det_ref_input(1,1,1),N_int,href)
+ i_pert_count = 0
+ do i = 1, n_det_outer_input
+  call i_h_j(psi_det_outer_input(1,1,i),psi_det_outer_input(1,1,i),N_int,hka)
+  f = 1.d0/(href - hka)
+  H_array = 0.d0
+  accu = 0.d0
+! n_h_out = number_of_holes(psi_det_outer_input(1,1,i))
+! n_p_out = number_of_particles(psi_det_outer_input(1,1,i))
+  do j=1,n_det_ref_input
+   n_h_in = number_of_holes(psi_det_ref_input(1,1,j))
+   n_p_in = number_of_particles(psi_det_ref_input(1,1,j))
+!  if(n_h_in == 0 .and. n_h_in == 0)then
+    call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka)
+!  else 
+!   hka = 0.d0
+!  endif
+   H_array(j) = hka
+   accu += psi_coef_ref_input(j,1) * hka
+  enddo
+  lambda_i = psi_coef_outer_input(i,1)/accu
+  i_pert = 1
+  if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then
+   i_pert = 0
+  endif
+  do j = 1, n_det_ref_input
+   if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
+    i_pert = 1
+    exit
+   endif
+  enddo
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! i_pert = 0
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  if(i_pert==1)then
+   lambda_i = f
+   i_pert_count +=1
+  endif
+  do k=1,n_det_ref_input
+    double precision :: contrib
+    contrib = H_array(k) * H_array(k) * lambda_i
+    dressing_matrix(k, k) += contrib
+    do j=k+1,n_det_ref_input
+      contrib = H_array(k) * H_array(j) * lambda_i
+      dressing_matrix(k, j) += contrib
+      dressing_matrix(j, k) += contrib
+    enddo 
+  enddo
+ enddo
+end
+
+
 subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, &
                                                            psi_det_outer_input,psi_coef_outer_input,n_det_outer_input)
  use bitmasks
@@ -112,16 +221,17 @@ subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi
    accu += psi_coef_ref_input(j,1) * hka
   enddo
   lambda_i = psi_coef_outer_input(i,1)/accu
-  i_pert = 1
+  i_pert = 0
   if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then
    i_pert = 0
   endif
   do j = 1, n_det_ref_input
-   if(dabs(H_array(j)*lambda_i).gt.0.3d0)then
+   if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
     i_pert = 1
     exit
    endif
   enddo
+! i_pert = 0
   if(i_pert==1)then
    lambda_i = f
    i_pert_count +=1
@@ -170,114 +280,379 @@ subroutine diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_
 
 end
 
-subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p)
+subroutine give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p)
  use bitmasks
  implicit none
- integer, intent(out) :: n_det_1h1p, n_det_2h1p
+ integer, intent(in)  :: i_hole
+ integer, intent(out) :: n_det_1h1p, n_det_2h1p,n_det_extra_1h_or_1p
  integer :: i
  integer :: n_det_ref_restart_tmp,n_det_1h
  integer ::  number_of_holes,n_h, number_of_particles,n_p
+ logical :: is_the_hole_in_det
  n_det_ref_restart_tmp = 0
  n_det_1h = 0
  n_det_1h1p = 0
  n_det_2h1p = 0
+ n_det_extra_1h_or_1p = 0
  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))
    if(n_h == 0 .and. n_p == 0)then
     n_det_ref_restart_tmp +=1
    else if (n_h ==1 .and. n_p==0)then
-    n_det_1h +=1
+    if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then  
+     n_det_1h +=1
+    else  
+     n_det_extra_1h_or_1p +=1
+    endif
+   else if (n_h ==0 .and. n_p==1)then
+    n_det_extra_1h_or_1p +=1
    else if (n_h ==1 .and. n_p==1)then
     n_det_1h1p +=1
    else if (n_h ==2 .and. n_p==1)then
     n_det_2h1p +=1
    else 
     print*,'PB !!!!'
-    print*,'You have something else than a 1h, 1h1p or 2h1p'
+    print*,'You have something else than a 1h, 1p, 1h1p or 2h1p'
+    print*,'n_h,n_p = ',n_h,n_p
     call debug_det(psi_det(1,1,i),N_int)
     stop
    endif
   enddo
-! if(n_det_1h.ne.1)then
-!  print*,'PB !! You have more than one 1h'
-!  stop
-! endif
  if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then
     print*,'PB !!!!'
     print*,'You have forgotten something in your generators ... '
     stop 
  endif
-
+ if(n_det_2h1p + n_det_1h1p + n_det_extra_1h_or_1p + n_det_generators .ne. N_det)then
+    print*,'PB !!!!'
+    print*,'You have forgotten something in your generators ... '
+    stop 
+ endif
 
 end
 
-subroutine give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p)
+subroutine give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p)
  use bitmasks
  implicit none
- integer, intent(out) :: n_det_1h1p, n_det_1h2p
+ integer, intent(out) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
  integer :: i
  integer :: n_det_ref_restart_tmp,n_det_1h
  integer ::  number_of_holes,n_h, number_of_particles,n_p
- n_det_ref_restart_tmp = 0
+ logical :: is_the_hole_in_det
- n_det_1h = 0
+ n_det_ref_1h_1p = 0
+ n_det_2h1p = 0
  n_det_1h1p = 0
- n_det_1h2p = 0
  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))
    if(n_h == 0 .and. n_p == 0)then
-    n_det_ref_restart_tmp +=1
+    n_det_ref_1h_1p +=1
+   else if (n_h ==1 .and. n_p==0)then
+     n_det_ref_1h_1p +=1
    else if (n_h ==0 .and. n_p==1)then
-    n_det_1h +=1
+    n_det_ref_1h_1p +=1
+   else if (n_h ==1 .and. n_p==1)then
+    n_det_1h1p +=1
+   else if (n_h ==2 .and. n_p==1)then
+    n_det_2h1p +=1
+   else 
+    print*,'PB !!!!'
+    print*,'You have something else than a 1h, 1p, 1h1p or 2h1p'
+    print*,'n_h,n_p = ',n_h,n_p
+    call debug_det(psi_det(1,1,i),N_int)
+    stop
+   endif
+  enddo
+
+end
+
+subroutine give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p)
+ use bitmasks
+ implicit none
+ integer, intent(out) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
+ integer :: i
+ integer :: n_det_ref_restart_tmp,n_det_1h
+ integer ::  number_of_holes,n_h, number_of_particles,n_p
+ logical :: is_the_hole_in_det
+ n_det_ref_1h_1p = 0
+ n_det_1h2p = 0
+ n_det_1h1p = 0
+ 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))
+   if(n_h == 0 .and. n_p == 0)then
+    n_det_ref_1h_1p +=1
+   else if (n_h ==1 .and. n_p==0)then
+     n_det_ref_1h_1p +=1
+   else if (n_h ==0 .and. n_p==1)then
+    n_det_ref_1h_1p +=1
    else if (n_h ==1 .and. n_p==1)then
     n_det_1h1p +=1
    else if (n_h ==1 .and. n_p==2)then
     n_det_1h2p +=1
    else 
     print*,'PB !!!!'
-    print*,'You have something else than a 1p, 1h1p or 1h2p'
+    print*,'You have something else than a 1h, 1p, 1h1p or 1h2p'
+    print*,'n_h,n_p = ',n_h,n_p
     call debug_det(psi_det(1,1,i),N_int)
     stop
    endif
   enddo
- if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then
+
+end
+
+subroutine give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& 
+                                                          psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p)
+ use bitmasks
+ implicit none
+ integer, intent(in) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
+ integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
+ integer(bit_kind), intent(out) :: psi_det_2h1p(N_int,2,n_det_2h1p)
+ integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
+ double precision, intent(out)  :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
+ double precision, intent(out)  :: psi_coef_2h1p(n_det_2h1p,N_states)
+ double precision, intent(out)  :: psi_coef_1h1p(n_det_1h1p,N_states)
+ integer :: n_det_ref_1h_1p_tmp,n_det_2h1p_tmp,n_det_1h1p_tmp
+ integer :: i,j
+ integer :: n_det_ref_restart_tmp,n_det_1h
+ integer ::  number_of_holes,n_h, number_of_particles,n_p
+ logical :: is_the_hole_in_det
+ integer, allocatable :: index_ref_1h_1p(:)
+ integer, allocatable :: index_2h1p(:)
+ integer, allocatable :: index_1h1p(:)
+ allocate(index_ref_1h_1p(n_det))
+ allocate(index_2h1p(n_det))
+ allocate(index_1h1p(n_det))
+ n_det_ref_1h_1p_tmp = 0
+ n_det_2h1p_tmp = 0
+ n_det_1h1p_tmp = 0
+ 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))
+   if(n_h == 0 .and. n_p == 0)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==1 .and. n_p==0)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==0 .and. n_p==1)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==1 .and. n_p==1)then
+    n_det_1h1p_tmp +=1
+    index_1h1p(n_det_1h1p_tmp) = i
+   else if (n_h ==2 .and. n_p==1)then
+    n_det_2h1p_tmp +=1
+    index_2h1p(n_det_2h1p_tmp) = i
+   else 
     print*,'PB !!!!'
-    print*,'You have forgotten something in your generators ... '
+    print*,'You have something else than a 1h, 1p, 1h1p or 2h1p'
-    stop 
+    print*,'n_h,n_p = ',n_h,n_p
- endif
+    call debug_det(psi_det(1,1,i),N_int)
+    stop
+   endif
+  enddo
+ do i = 1, n_det_2h1p
+  do j = 1, N_int
+   psi_det_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i))
+   psi_det_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j)
+  enddo
+ enddo
+
+ do i = 1, n_det_1h1p
+  do j = 1, N_int
+   psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
+   psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
+  enddo
+ enddo
+
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, N_int
+   psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
+   psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
+  enddo
+ enddo
+
+end
+
+subroutine give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& 
+                                                          psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p)
+ use bitmasks
+ implicit none
+ integer, intent(in) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
+ integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
+ integer(bit_kind), intent(out) :: psi_det_1h2p(N_int,2,n_det_1h2p)
+ integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
+ double precision, intent(out)  :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
+ double precision, intent(out)  :: psi_coef_1h2p(n_det_1h2p,N_states)
+ double precision, intent(out)  :: psi_coef_1h1p(n_det_1h1p,N_states)
+ integer :: n_det_ref_1h_1p_tmp,n_det_1h2p_tmp,n_det_1h1p_tmp
+ integer :: i,j
+ integer :: n_det_ref_restart_tmp,n_det_1h
+ integer ::  number_of_holes,n_h, number_of_particles,n_p
+ logical :: is_the_hole_in_det
+ integer, allocatable :: index_ref_1h_1p(:)
+ integer, allocatable :: index_1h2p(:)
+ integer, allocatable :: index_1h1p(:)
+ allocate(index_ref_1h_1p(n_det))
+ allocate(index_1h2p(n_det))
+ allocate(index_1h1p(n_det))
+ n_det_ref_1h_1p_tmp = 0
+ n_det_1h2p_tmp = 0
+ n_det_1h1p_tmp = 0
+ 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))
+   if(n_h == 0 .and. n_p == 0)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==1 .and. n_p==0)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==0 .and. n_p==1)then
+    n_det_ref_1h_1p_tmp +=1
+    index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
+   else if (n_h ==1 .and. n_p==1)then
+    n_det_1h1p_tmp +=1
+    index_1h1p(n_det_1h1p_tmp) = i
+   else if (n_h ==1 .and. n_p==2)then
+    n_det_1h2p_tmp +=1
+    index_1h2p(n_det_1h2p_tmp) = i
+   else 
+    print*,'PB !!!!'
+    print*,'You have something else than a 1h, 1p, 1h1p or 1h2p'
+    print*,'n_h,n_p = ',n_h,n_p
+    call debug_det(psi_det(1,1,i),N_int)
+    stop
+   endif
+  enddo
+ do i = 1, n_det_1h2p
+  do j = 1, N_int
+   psi_det_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i))
+   psi_det_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j)
+  enddo
+ enddo
+
+ do i = 1, n_det_1h1p
+  do j = 1, N_int
+   psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
+   psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
+  enddo
+ enddo
+
+ do i = 1, n_det_ref_1h_1p
+  do j = 1, N_int
+   psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
+   psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
+  enddo
+ enddo
+
+end
+
+
+
+subroutine give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p)
+ use bitmasks
+ implicit none
+ integer, intent(in)  ::i_particl
+ integer, intent(out) :: n_det_1h1p, n_det_1h2p,n_det_extra_1h_or_1p
+ integer :: i
+ integer :: n_det_ref_restart_tmp,n_det_1p
+ integer ::  number_of_holes,n_h, number_of_particles,n_p
+ logical :: is_the_particl_in_det
+ n_det_ref_restart_tmp = 0
+ n_det_1p = 0
+ n_det_1h1p = 0
+ n_det_1h2p = 0
+ n_det_extra_1h_or_1p = 0
+ 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))
+   if(n_h == 0 .and. n_p == 0)then
+    n_det_ref_restart_tmp +=1
+   else if (n_h ==0 .and. n_p==1)then
+    if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then
+     n_det_1p +=1
+    else 
+     n_det_extra_1h_or_1p +=1
+    endif
+   else if (n_h ==1 .and. n_p==0)then
+    n_det_extra_1h_or_1p +=1
+   else if (n_h ==1 .and. n_p==1)then
+    n_det_1h1p +=1
+   else if (n_h ==1 .and. n_p==2)then
+    n_det_1h2p +=1
+   else 
+    print*,'PB !!!!'
+    print*,'You have something else than a 1h, 1p, 1h1p or 1h2p'
+    call debug_det(psi_det(1,1,i),N_int)
+    stop
+   endif
+  enddo
+!if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then
+!   print*,'PB !!!!'
+!   print*,'You have forgotten something in your generators ... '
+!   stop 
+!endif
 
 
 end
 
 
-subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p)
+subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
  use bitmasks
  implicit none
- integer, intent(in) :: n_det_1h1p,n_det_2h1p
+ integer, intent(in) :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p,i_hole
  integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators)
  integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p)
  integer(bit_kind), intent(out) :: psi_2h1p(N_int,2,n_det_2h1p)
+ integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)
  double precision,  intent(out) :: psi_ref_coef_out(N_det_generators,N_states)
  double precision,  intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states)
  double precision,  intent(out) :: psi_coef_2h1p(n_det_2h1p, N_states)
+ double precision,  intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states)
 
  integer :: i,j
  integer :: degree
  integer ::  number_of_holes,n_h, number_of_particles,n_p
- integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp
+ integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp,n_det_extra_1h_or_1p_tmp
+ integer :: n_det_1h_tmp
  integer, allocatable  :: index_generator(:)
  integer, allocatable :: index_1h1p(:)
  integer, allocatable :: index_2h1p(:)
+ integer, allocatable :: index_extra_1h_or_1p(:)
+ logical :: is_the_hole_in_det
 
  allocate(index_1h1p(n_det))
  allocate(index_2h1p(n_det))
+ allocate(index_extra_1h_or_1p(n_det))
  allocate(index_generator(N_det))
 
 
  n_det_generators_tmp = 0
  n_det_1h1p_tmp = 0
  n_det_2h1p_tmp = 0
+ n_det_extra_1h_or_1p_tmp = 0
+ n_det_1h_tmp = 0
  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))
@@ -287,6 +662,16 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o
    else if (n_h ==2 .and. n_p==1)then
     n_det_2h1p_tmp +=1
     index_2h1p(n_det_2h1p_tmp) = i
+   else if (n_h ==0 .and. n_p==1)then
+    n_det_extra_1h_or_1p_tmp +=1
+    index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i
+   else if (n_h ==1 .and. n_p==0)then
+    if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then  
+     n_det_1h_tmp +=1
+    else
+     n_det_extra_1h_or_1p_tmp +=1
+     index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i
+    endif
    endif
    do j = 1, N_det_generators
     call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int)
@@ -315,6 +700,12 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o
   stop
  endif
 
+ if(n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp)then
+  print*,'PB !!!'
+  print*,'n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp'
+  stop
+ endif
+
  do i = 1,N_det_generators
   do j = 1, N_int
    psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i))
@@ -345,41 +736,59 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o
   enddo
  enddo
 
+ do i = 1, n_det_extra_1h_or_1p
+  do j = 1, N_int
+   psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i))
+   psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j)
+  enddo
+ enddo
 
  deallocate(index_generator)
  deallocate(index_1h1p)
  deallocate(index_2h1p)
+ deallocate(index_extra_1h_or_1p)
 
 end
 
 
-subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p)
+subroutine split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
  use bitmasks
  implicit none
- integer, intent(in) :: n_det_1h1p,n_det_1h2p
+ integer, intent(in) :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p,i_particl
  integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators)
  integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p)
  integer(bit_kind), intent(out) :: psi_1h2p(N_int,2,n_det_1h2p)
+ integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)
  double precision,  intent(out) :: psi_ref_coef_out(N_det_generators,N_states)
  double precision,  intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states)
  double precision,  intent(out) :: psi_coef_1h2p(n_det_1h2p, N_states)
+ double precision,  intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states)
 
  integer :: i,j
  integer :: degree
  integer ::  number_of_holes,n_h, number_of_particles,n_p
- integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp
+ integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp,n_det_extra_1h_or_1p_tmp
  integer, allocatable  :: index_generator(:)
  integer, allocatable :: index_1h1p(:)
  integer, allocatable :: index_1h2p(:)
+ integer, allocatable :: index_extra_1h_or_1p(:)
+ logical :: is_the_particl_in_det
+ integer :: n_det_1p_tmp
 
  allocate(index_1h1p(n_det))
  allocate(index_1h2p(n_det))
+ allocate(index_extra_1h_or_1p(n_det))
  allocate(index_generator(N_det))
 
 
  n_det_generators_tmp = 0
  n_det_1h1p_tmp = 0
  n_det_1h2p_tmp = 0
+ n_det_extra_1h_or_1p_tmp = 0
+ n_det_1p_tmp = 0
  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))
@@ -389,6 +798,15 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o
    else if (n_h ==1 .and. n_p==2)then
     n_det_1h2p_tmp +=1
     index_1h2p(n_det_1h2p_tmp) = i
+   else if (n_h ==1 .and. n_p==0)then
+    n_det_extra_1h_or_1p_tmp +=1
+    index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i
+   else if (n_h ==0 .and. n_p==1)then
+    if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then
+     n_det_1p_tmp +=1
+    else 
+     n_det_extra_1h_or_1p_tmp +=1
+    endif
    endif
    do j = 1, N_det_generators
     call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int)
@@ -448,9 +866,20 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o
  enddo
 
 
+ do i = 1, n_det_extra_1h_or_1p
+  do j = 1, N_int
+   psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i))
+   psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i))
+  enddo
+  do j = 1, N_states
+   psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j)
+  enddo
+ enddo
+
  deallocate(index_generator)
  deallocate(index_1h1p)
  deallocate(index_1h2p)
+ deallocate(index_extra_1h_or_1p)
 
 end
 

plugins/FOBOCI/routines_foboci.irp.f

@@ -332,20 +332,20 @@ subroutine save_osoci_natural_mos
  enddo
 
  tmp = tmp_bis
-!! Symetrization act-virt
+!!! Symetrization act-virt
- do j = 1, n_virt_orb
+! do j = 1, n_virt_orb
-  j_virt= list_virt(j)
+!  j_virt= list_virt(j)
-  accu = 0.d0
+!  accu = 0.d0
-  do i = 1, n_act_orb
+!  do i = 1, n_act_orb
-   jorb = list_act(i)
+!   jorb = list_act(i)
-   accu += dabs(tmp_bis(j_virt,jorb))
+!   accu += dabs(tmp_bis(j_virt,jorb))
-  enddo
+!  enddo
-  do i = 1, n_act_orb
+!  do i = 1, n_act_orb
-   iorb = list_act(i)
+!   iorb = list_act(i)
-   tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
+!   tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
-   tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
+!   tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
-  enddo
+!  enddo
- enddo
+! enddo
 
 !! Symetrization act-inact
 !do j = 1, n_inact_orb
@@ -387,16 +387,16 @@ 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,dabs(tmp(iorb,jorb))
+    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,dabs(tmp(iorb,jorb))
+    print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
   enddo
  enddo
@@ -410,8 +410,9 @@ subroutine save_osoci_natural_mos
  enddo
 
  label = "Natural"
+ 
  call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
- soft_touch mo_coef
+!soft_touch mo_coef
  deallocate(tmp,occ)
 
 
@@ -518,16 +519,16 @@ 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,dabs(tmp(iorb,jorb))
+    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,dabs(tmp(iorb,jorb))
+    print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
    endif
   enddo
  enddo
@@ -602,15 +603,210 @@ end
  
  subroutine provide_properties
    implicit none
-   integer :: i
+   call print_mulliken_sd
-   double precision :: accu
+   call print_hcc
-   if(.True.)then
-    accu= 0.d0
-    do i = 1, nucl_num
-     accu += mulliken_spin_densities(i)
-     print*,i,nucl_charge(i),mulliken_spin_densities(i)
-    enddo
-    print*,'Sum of Mulliken SD = ',accu
-   endif
  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

plugins/Generators_restart/generators.irp.f

@@ -1,5 +1,5 @@
 use bitmasks
-
+ 
 BEGIN_PROVIDER [ integer, N_det_generators ]
  implicit none
  BEGIN_DOC
@@ -8,17 +8,18 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  integer :: i
  integer, save :: ifirst = 0
  double precision :: norm
- read_wf = .True.
  if(ifirst == 0)then
-  N_det_generators = N_det
+  call ezfio_get_determinants_n_det(N_det_generators)
   ifirst = 1
+ else
+  print*,'PB in generators restart !!!'
  endif
  call write_int(output_determinants,N_det_generators,'Number of generators')
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ]
-&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ]
  implicit none
  BEGIN_DOC
  ! read wf
@@ -26,17 +27,20 @@ END_PROVIDER
  END_DOC
  integer                        :: i, k
  integer, save :: ifirst = 0
+ double precision, allocatable  :: psi_coef_read(:,:)
  if(ifirst == 0)then
-  do i=1,N_det_generators
+  call read_dets(psi_det_generators,N_int,N_det_generators)
-    do k=1,N_int
+   allocate (psi_coef_read(N_det_generators,N_states))
-      psi_det_generators(k,1,i) = psi_det(k,1,i)
+   call ezfio_get_determinants_psi_coef(psi_coef_read)
-      psi_det_generators(k,2,i) = psi_det(k,2,i)
-    enddo
    do k = 1, N_states
-    psi_coef_generators(i,k) = psi_coef(i,k)
+    do i = 1, N_det_generators
+     psi_coef_generators(i,k) = psi_coef_read(i,k)
+    enddo
    enddo
-  enddo
   ifirst = 1
+  deallocate(psi_coef_read)
+ else 
+  print*,'PB in generators restart !!!'
  endif
 
 END_PROVIDER

plugins/Hartree_Fock/damping_SCF.irp.f

@@ -119,7 +119,9 @@ subroutine damping_SCF
   write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  '====','================','================','================', '===='
   write(output_hartree_fock,*)
   
-  call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+  if(.not.no_oa_or_av_opt)then
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+  endif
 
   call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
   call ezfio_set_hartree_fock_energy(E_min)

plugins/Perturbation/pt2_equations.irp.f

@@ -126,6 +126,8 @@ subroutine pt2_moller_plesset ($arguments)
     delta_e = (Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)) + &
               (Fock_matrix_diag_mo(h2) - Fock_matrix_diag_mo(p2))
     delta_e = 1.d0/delta_e
+!   print*,'h1,p1',h1,p1
+!   print*,'h2,p2',h2,p2
   else if (degree == 1) then
     call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
     delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1) 

plugins/Properties/hyperfine_constants.irp.f

@@ -133,3 +133,16 @@ END_PROVIDER
  enddo
 
 END_PROVIDER 
+
+
+subroutine print_hcc
+ implicit none
+ double precision :: accu
+ integer :: i,j
+ print*,'Z               AU           GAUSS              MHZ             cm^-1'
+ do i = 1, nucl_num
+  write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
+ enddo
+
+end
+

plugins/Properties/mulliken.irp.f

@@ -105,3 +105,34 @@ END_PROVIDER
  enddo
 
 END_PROVIDER
+
+
+subroutine print_mulliken_sd
+ implicit none
+ double precision :: accu
+ integer :: i
+ integer :: j
+ print*,'Mulliken spin densities'
+ accu= 0.d0
+ do i = 1, nucl_num
+  print*,i,nucl_charge(i),mulliken_spin_densities(i)
+  accu += mulliken_spin_densities(i)
+ enddo
+ print*,'Sum of Mulliken SD = ',accu
+ print*,'AO SPIN POPULATIONS'
+ accu = 0.d0
+ do i = 1, ao_num
+  accu += spin_gross_orbital_product(i)
+  write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
+ enddo
+ print*,'sum = ',accu
+ accu = 0.d0
+ print*,'Angular momentum analysis'
+ do i = 0,  ao_l_max
+  accu += spin_population_angular_momentum(i)
+  print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
+ print*,'sum = ',accu
+ enddo
+
+end
+

plugins/Properties/print_hcc.irp.f

@@ -1,17 +1,6 @@
-program print_hcc
+program print_hcc_main
  implicit none
  read_wf = .True.
  touch read_wf
- call test
+ call print_hcc
 end
-subroutine test
- implicit none
- double precision :: accu
- integer :: i,j
- print*,'Z               AU           GAUSS              MHZ             cm^-1'
- do i = 1, nucl_num
-  write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
- enddo
-
-end
-

plugins/Properties/print_mulliken.irp.f

@@ -2,34 +2,5 @@ program print_mulliken
  implicit none
  read_wf = .True.
  touch read_wf
- print*,'Mulliken spin densities'
+ call print_mulliken_sd
-
- call test
 end
-subroutine test
- double precision :: accu
- integer :: i
- integer :: j
- accu= 0.d0
- do i = 1, nucl_num
-  print*,i,nucl_charge(i),mulliken_spin_densities(i)
-  accu += mulliken_spin_densities(i)
- enddo
- print*,'Sum of Mulliken SD = ',accu
- print*,'AO SPIN POPULATIONS'
- accu = 0.d0
- do i = 1, ao_num
-  accu += spin_gross_orbital_product(i)
-  write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
- enddo
- print*,'sum = ',accu
- accu = 0.d0
- print*,'Angular momentum analysis'
- do i = 0,  ao_l_max
-  accu += spin_population_angular_momentum(i)
-  print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
- print*,'sum = ',accu
- enddo
-
-end
-

scripts/generate_h_apply.py

@@ -8,11 +8,22 @@ copy_buffer
 declarations
 decls_main
 deinit_thread
-do_double_excitations
+skip
+init_main
+filter_integrals
+filter2p
+filter2h2p_double
+filter2h2p_single
 filter1h
 filter1p
-filter2h2p
+only_2p_single
-filter2p
+only_2p_double
+filter_only_1h1p_single
+filter_only_1h1p_double
+filter_only_1h2p_single
+filter_only_1h2p_double
+filter_only_2h2p_single
+filter_only_2h2p_double
 filterhole
 filter_integrals
 filter_only_1h1p_double
@@ -181,7 +192,7 @@ class H_apply(object):
      if (is_a_2p(hole)) cycle
     """
   def filter_1p(self):
-    self["filter0p"] = """
+    self["filter1p"] = """
 !    ! DIR$ FORCEINLINE
      if (is_a_1p(hole)) cycle
     """
@@ -207,6 +218,27 @@ class H_apply(object):
      if (is_a_1h1p(key).eqv..False.) cycle
     """
 
+  def filter_only_2h2p(self):
+    self["filter_only_2h2p_single"] = """
+!    ! DIR$ FORCEINLINE
+     if (is_a_two_holes_two_particles(hole).eqv..False.) cycle
+    """
+    self["filter_only_1h1p_double"] = """
+!    ! DIR$ FORCEINLINE
+     if (is_a_two_holes_two_particles(key).eqv..False.) cycle
+    """
+
+
+  def filter_only_1h2p(self):
+    self["filter_only_1h2p_single"] = """
+!    ! DIR$ FORCEINLINE
+     if (is_a_1h2p(hole).eqv..False.) cycle
+    """
+    self["filter_only_1h2p_double"] = """
+!    ! DIR$ FORCEINLINE
+     if (is_a_1h2p(key).eqv..False.) cycle
+    """
+
 
   def unset_skip(self):
     self["skip"] = """
@@ -214,9 +246,12 @@ class H_apply(object):
 
 
   def set_filter_2h_2p(self):
-    self["filter2h2p"] = """
+    self["filter2h2p_double"] = """
      if (is_a_two_holes_two_particles(key)) cycle
     """
+    self["filter2h2p_single"] = """
+     if (is_a_two_holes_two_particles(hole)) cycle
+    """
 
 
   def set_perturbation(self,pert):

src/Bitmask/bitmask_cas_routines.irp.f

@@ -212,6 +212,12 @@ logical function is_a_two_holes_two_particles(key_in)
  implicit none
  integer(bit_kind), intent(in) :: key_in(N_int,2)
  integer :: i,i_diff
+ integer :: number_of_holes, number_of_particles
+ is_a_two_holes_two_particles = .False.
+ if(number_of_holes(key_in) == 2 .and. number_of_particles(key_in) == 2)then
+  is_a_two_holes_two_particles = .True. 
+  return
+ endif
  i_diff = 0
  if(N_int == 1)then
    i_diff = i_diff  &
@@ -456,6 +462,17 @@ logical function is_a_1h1p(key_in)
 
 end
 
+logical function is_a_1h2p(key_in)
+ implicit none
+ integer(bit_kind), intent(in) :: key_in(N_int,2)
+ integer :: number_of_particles, number_of_holes
+ is_a_1h2p = .False.
+ if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.2)then
+  is_a_1h2p = .True.
+ endif
+
+end
+
 logical function is_a_1h(key_in)
  implicit none
  integer(bit_kind), intent(in) :: key_in(N_int,2)

src/Bitmask/bitmasks.irp.f

@@ -95,9 +95,40 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ]
+ implicit none
+ BEGIN_DOC
+ ! Number of bitmasks for generators
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ 
+ call ezfio_has_bitmasks_N_mask_gen(exists)
+ if (exists) then
+   call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart)
+   integer                        :: N_int_check
+   integer                        :: bit_kind_check
+   call ezfio_get_bitmasks_bit_kind(bit_kind_check)
+   if (bit_kind_check /= bit_kind) then
+     print *,  bit_kind_check, bit_kind
+     print *,  'Error: bit_kind is not correct in EZFIO file'
+   endif
+   call ezfio_get_bitmasks_N_int(N_int_check)
+   if (N_int_check /= N_int) then
+     print *,  N_int_check, N_int
+     print *,  'Error: N_int is not correct in EZFIO file'
+   endif
+ else
+   N_generators_bitmask_restart = 1
+ endif
+ ASSERT (N_generators_bitmask_restart > 0)
+
+END_PROVIDER
 
 
-BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask) ]
+
+
+BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ]
  implicit none
  BEGIN_DOC
  ! Bitmasks for generator determinants.
@@ -306,7 +337,7 @@ END_PROVIDER
 
  n_inact_orb = 0
  n_virt_orb = 0
- if(N_generators_bitmask == 1)then
+ if(N_generators_bitmask_restart == 1)then
   do j = 1, N_int
    inact_bitmask(j,1) = xor(generators_bitmask_restart(j,1,1,1),cas_bitmask(j,1,1))
    inact_bitmask(j,2) = xor(generators_bitmask_restart(j,2,1,1),cas_bitmask(j,2,1))
@@ -319,15 +350,15 @@ END_PROVIDER
    i_hole = 1
    i_gen = 1
    do i = 1, N_int
-     inact_bitmask(i,1) = generators_bitmask(i,1,i_hole,i_gen)
+     inact_bitmask(i,1) = generators_bitmask_restart(i,1,i_hole,i_gen)
-     inact_bitmask(i,2) = generators_bitmask(i,2,i_hole,i_gen)
+     inact_bitmask(i,2) = generators_bitmask_restart(i,2,i_hole,i_gen)
      n_inact_orb += popcnt(inact_bitmask(i,1))
    enddo
    i_part = 2
    i_gen = 3
    do i = 1, N_int
-     virt_bitmask(i,1) = generators_bitmask(i,1,i_part,i_gen)
+     virt_bitmask(i,1) = generators_bitmask_restart(i,1,i_part,i_gen)
-     virt_bitmask(i,2) = generators_bitmask(i,2,i_part,i_gen)
+     virt_bitmask(i,2) = generators_bitmask_restart(i,2,i_part,i_gen)
      n_virt_orb  += popcnt(virt_bitmask(i,1))
    enddo
  endif

src/Determinants/H_apply.template.f

@@ -175,6 +175,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
   
   logical :: check_double_excitation 
   logical :: is_a_1h1p
+  logical :: is_a_1h2p
   logical :: is_a_1h
   logical :: is_a_1p
   logical :: is_a_2p
@@ -304,8 +305,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
           k = ishft(j_b-1,-bit_kind_shift)+1
           l = j_b-ishft(k-1,bit_kind_shift)-1
           key(k,other_spin) = ibset(key(k,other_spin),l)
-          $filter2h2p
+          $filter2h2p_double
           $filter_only_1h1p_double
+          $filter_only_1h2p_double
+          $filter_only_2h2p_double
           $only_2p_double
           key_idx += 1
           do k=1,N_int
@@ -353,8 +356,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
         k = ishft(j_b-1,-bit_kind_shift)+1
         l = j_b-ishft(k-1,bit_kind_shift)-1
         key(k,ispin) = ibset(key(k,ispin),l)
-        $filter2h2p
+        $filter2h2p_double
         $filter_only_1h1p_double
+        $filter_only_1h2p_double
+        $filter_only_2h2p_double
         $only_2p_double
         key_idx += 1
         do k=1,N_int
@@ -424,6 +429,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   
   logical :: check_double_excitation 
   logical :: is_a_1h1p
+  logical :: is_a_1h2p
   logical :: is_a_1h
   logical :: is_a_1p
   logical :: is_a_2p
@@ -505,8 +511,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
       $filter1h
       $filter1p
       $filter2p
-      $filter2h2p
+      $filter2h2p_single
       $filter_only_1h1p_single
+      $filter_only_1h2p_single
+      $filter_only_2h2p_single
       key_idx += 1
       do k=1,N_int
         keys_out(k,1,key_idx) = hole(k,1)
@@ -532,4 +540,3 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   
 end
 
-

src/Determinants/SC2.irp.f

@@ -1,4 +1,4 @@
-subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
+subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -21,6 +21,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
   integer(bit_kind), intent(in)  :: dets_in(Nint,2,sze)
   double precision, intent(inout) :: u_in(dim_in,N_st)
   double precision, intent(out)  :: energies(N_st)
+  double precision, intent(out)  :: diag_H_elements(dim_in)
   double precision, intent(in)   :: convergence
   ASSERT (N_st > 0)
   ASSERT (sze > 0)
@@ -197,6 +198,9 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
     converged =  dabs(e_corr_double - e_corr_double_before) < convergence
     converged = converged 
     if (converged) then
+      do i = 1, dim_in
+       diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i)
+      enddo
       exit
     endif
     e_corr_double_before = e_corr_double

src/Determinants/determinants.irp.f

@@ -58,7 +58,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
  else
    psi_det_size = 1
  endif
- psi_det_size = max(psi_det_size,10000)
+ psi_det_size = max(psi_det_size,100000)
  call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays')
 
 END_PROVIDER

src/Determinants/diagonalize_CI_SC2.irp.f

@@ -23,8 +23,10 @@ END_PROVIDER
  threshold_convergence_SC2 = 1.d-10
 
  END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
 &BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ]
   implicit none
   BEGIN_DOC
   ! Eigenvectors/values of the CI matrix
@@ -39,7 +41,8 @@ END_PROVIDER
   enddo
   
   call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
-        size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+!       size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+        diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
 END_PROVIDER
 
 subroutine diagonalize_CI_SC2
@@ -54,5 +57,6 @@ subroutine diagonalize_CI_SC2
       psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
     enddo
   enddo
-  SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors
+  SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2
+! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors 
 end

src/Determinants/save_natorb.irp.f

@@ -2,5 +2,6 @@ program save_natorb
   read_wf = .True.
   touch read_wf
   call save_natural_mos
+  call save_ref_determinant
 end
 

src/Integrals_Bielec/map_integrals.irp.f

@@ -230,7 +230,6 @@ subroutine clear_ao_map
 end
 
 
-
 !! MO Map
 !! ======
 

src/Integrals_Bielec/mo_bi_integrals.irp.f

@@ -72,7 +72,7 @@ subroutine add_integrals_to_map(mask_ijkl)
   integer                        :: i2,i3,i4
   double precision,parameter     :: thr_coef = 1.d-10
   
-  PROVIDE ao_bielec_integrals_in_map 
+  PROVIDE ao_bielec_integrals_in_map  mo_coef
   
   !Get list of MOs for i,j,k and l
   !-------------------------------
@@ -329,7 +329,7 @@ end
   double precision, allocatable  :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
   
   if (.not.do_direct_integrals) then
-    PROVIDE ao_bielec_integrals_in_map
+    PROVIDE ao_bielec_integrals_in_map mo_coef
   endif
   
   mo_bielec_integral_jj_from_ao = 0.d0
@@ -495,4 +495,13 @@ subroutine clear_mo_map
   call map_deinit(mo_integrals_map)
   FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
   FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
+
+
+end
+
+subroutine provide_all_mo_integrals
+ implicit none
+ provide mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
+ provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
+
 end

src/Integrals_Monoelec/mo_mono_ints.irp.f

@@ -5,6 +5,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
  ! array of the mono electronic hamiltonian on the MOs basis
  ! : sum of the kinetic and nuclear electronic potential 
   END_DOC
+  print*,'Providing the mono electronic integrals'
   do j = 1, mo_tot_num
    do i = 1, mo_tot_num
     mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)