diff --git a/config/ifort.cfg b/config/ifort.cfg
index f8c36e4e..cc848cba 100644
--- a/config/ifort.cfg
+++ b/config/ifort.cfg
@@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
 # 0 : Deactivate
 # 
 [OPTION]
-MODE    : DEBUG      ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
 CACHE   : 1          ; Enable cache_compile.py
 OPENMP  : 1          ; Append OpenMP flags
 
diff --git a/plugins/Hartree_Fock/.gitignore b/plugins/Hartree_Fock/.gitignore
index 9f1c0929..f1a4ff4f 100644
--- a/plugins/Hartree_Fock/.gitignore
+++ b/plugins/Hartree_Fock/.gitignore
@@ -5,7 +5,6 @@ AO_Basis
 Bitmask
 Electrons
 Ezfio_files
-Huckel_guess
 IRPF90_man
 IRPF90_temp
 Integrals_Bielec
@@ -16,7 +15,6 @@ Makefile
 Makefile.depend
 Nuclei
 Pseudo
-SCF
 Utils
 ZMQ
 ezfio_interface.irp.f
diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg
index d8207cc4..2fa29cf0 100644
--- a/plugins/Hartree_Fock/EZFIO.cfg
+++ b/plugins/Hartree_Fock/EZFIO.cfg
@@ -26,3 +26,10 @@ default: Huckel
 type: double precision
 doc: Calculated HF energy
 interface: ezfio
+
+[no_oa_or_av_opt]
+type: logical
+doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the SCF procedure
+interface: ezfio,provider,ocaml
+default: False
+
diff --git a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
index 85bdd3ad..6fb87e35 100644
--- a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
+++ b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Integrals_Bielec MOGuess 
+Integrals_Bielec MOGuess Bitmask
diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f
index 7c424aeb..c80077b3 100644
--- a/plugins/Hartree_Fock/diagonalize_fock.irp.f
+++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f
@@ -11,63 +11,35 @@
    double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
 
-!   if (mo_tot_num == ao_num) then
-!       ! Solve H.C = E.S.C in AO basis set
-!
-!       allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
-!       do j=1,ao_num
-!         do i=1,ao_num
-!           S(i,j) = ao_overlap(i,j)
-!           F(i,j) = Fock_matrix_ao(i,j)
-!         enddo
-!       enddo
-!
-!       n = ao_num
-!       lwork = 1+6*n + 2*n*n
-!       liwork = 3 + 5*n
-!       
-!       allocate(work(lwork), iwork(liwork) )
-!
-!       lwork = -1
-!       liwork = -1
-!
-!       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
-!         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
-!
-!       if (info /= 0) then
-!         print *,  irp_here//' failed : ', info
-!         stop 1
-!       endif
-!       lwork = int(work(1))
-!       liwork = iwork(1)
-!       deallocate(work,iwork)
-!       allocate(work(lwork), iwork(liwork) )
-!
-!       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
-!         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
-!
-!       if (info /= 0) then
-!         print *,  irp_here//' failed : ', info
-!         stop 1
-!       endif
-!       do j=1,mo_tot_num
-!         do i=1,ao_num
-!           eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
-!         enddo
-!       enddo
-!
-!       deallocate(work, iwork, F, S)
-!
-!  else 
-!
-       ! Solve H.C = E.C in MO basis set
-       
        allocate( F(mo_tot_num_align,mo_tot_num) )
        do j=1,mo_tot_num
          do i=1,mo_tot_num
            F(i,j) = Fock_matrix_mo(i,j)
          enddo
        enddo
+       if(no_oa_or_av_opt)then
+        integer :: iorb,jorb
+        do i = 1, n_act_orb
+         iorb = list_act(i)
+         do j = 1, n_inact_orb
+          jorb = list_inact(j)
+          F(iorb,jorb) = 0.d0
+          F(jorb,iorb) = 0.d0
+         enddo
+         do j = 1, n_virt_orb
+          jorb = list_virt(j)
+          F(iorb,jorb) = 0.d0
+          F(jorb,iorb) = 0.d0
+         enddo
+         do j = 1, n_core_orb
+          jorb = list_core(j)
+          F(iorb,jorb) = 0.d0
+          F(jorb,iorb) = 0.d0
+         enddo
+        enddo
+       endif
+
+  
        
 
        ! Insert level shift here
diff --git a/plugins/OVB/NEEDED_CHILDREN_MODULES b/plugins/OVB/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..05706cac
--- /dev/null
+++ b/plugins/OVB/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Determinants  Psiref_CAS
diff --git a/plugins/OVB/README.rst b/plugins/OVB/README.rst
new file mode 100644
index 00000000..13488c7e
--- /dev/null
+++ b/plugins/OVB/README.rst
@@ -0,0 +1,20 @@
+=======================
+OVB
+=======================
+The present module proposes an orthogonal Valence Bond analysis 
+of the wave function, that are the printing of the various Hamiltonian 
+matrix elements on the basis of the level of ionicity of the components 
+of the wave function. 
+
+Assumptions : it supposes that you have some orthogonal local orbitals within 
+the active space and that you performed a CI within the active orbitals. 
+Such CI might be complete or not, no matter. 
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
diff --git a/plugins/OVB/ovb_components.irp.f b/plugins/OVB/ovb_components.irp.f
new file mode 100644
index 00000000..82f0f931
--- /dev/null
+++ b/plugins/OVB/ovb_components.irp.f
@@ -0,0 +1,510 @@
+
+ use bitmasks
+ BEGIN_PROVIDER [integer, max_number_ionic]
+&BEGIN_PROVIDER [integer, min_number_ionic]
+ BEGIN_DOC
+ ! Maximum and minimum number of ionization in psi_ref 
+ END_DOC
+ implicit none
+ integer :: i,j
+ integer :: n_closed_shell_cas
+ max_number_ionic = 0
+ min_number_ionic = 100000
+ do i = 1, N_det_ref
+  j = n_closed_shell_cas(psi_ref(1,1,i),n_int)
+  if(j> max_number_ionic)then
+   max_number_ionic = j
+  endif
+  if(j< min_number_ionic)then
+   min_number_ionic = j
+  endif
+  
+ enddo
+ print*,'max_number_ionic = ',max_number_ionic
+ print*,'min_number_ionic = ',min_number_ionic
+END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, ionic_index, (min_number_ionic:max_number_ionic,0:N_det_ref)]
+&BEGIN_PROVIDER [double precision, normalization_factor_ionic, (min_number_ionic:max_number_ionic, N_states)]
+ BEGIN_DOC
+ ! Index of the various determinants in psi_ref according to their level of ionicity
+ ! ionic_index(i,0) = number of determinants in psi_ref having the degree of ionicity "i"
+ ! ionic_index(i,j) = index of the determinants having the degree of ionicity "i" 
+ END_DOC
+ implicit none
+ integer :: i,j,k
+ integer :: n_closed_shell_cas
+ double precision :: accu
+ ionic_index = 0
+ do i = 1, N_det_ref
+  j = n_closed_shell_cas(psi_ref(1,1,i),n_int)
+  ionic_index(j,0) +=1
+  ionic_index(j,ionic_index(j,0)) = i
+ enddo
+ do i = min_number_ionic,max_number_ionic 
+  accu = 0.d0
+  do j = 1, N_states
+   do k = 1, ionic_index(i,0)
+    accu += psi_ref_coef_diagonalized(ionic_index(i,k),j)   *  psi_ref_coef_diagonalized(ionic_index(i,k),j)
+   enddo
+   normalization_factor_ionic(i,j) = 1.d0/dsqrt(accu)
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [double precision, H_OVB_naked, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)]
+ BEGIN_DOC 
+ ! Hamiltonian matrix expressed in the basis of contracted forms in terms of ionic structures
+ END_DOC
+ implicit none
+ integer :: i,j,istate,k,l
+ double precision :: accu,hij
+ do i = min_number_ionic,max_number_ionic
+  do j = min_number_ionic,max_number_ionic
+   do istate = 1, N_states
+    accu = 0.d0
+    do k = 1,  ionic_index(i,0)
+     do l = 1,  ionic_index(j,0)
+      hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
+      accu += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * & 
+              psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
+     enddo
+    enddo
+    H_OVB_naked(i,j,istate) = accu
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [integer, n_couples_act_orb]
+ implicit none
+ n_couples_act_orb = 3
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, couples_act_orb, (n_couples_act_orb,2) ]
+ implicit none
+ 
+ couples_act_orb(1,1) = 20
+ couples_act_orb(1,2) = 21
+ couples_act_orb(2,1) = 22
+ couples_act_orb(2,2) = 23
+ couples_act_orb(3,1) = 24
+ couples_act_orb(3,2) = 25
+
+ END_PROVIDER 
+
+
+ BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom , (n_act_orb,n_act_orb)]
+ implicit none
+ BEGIN_DOC
+! Hamiltonian matrix elements between the various contracted functions 
+! that have a negative charge on a given active orbital
+ END_DOC
+ integer :: i,j,k,l,jj,ii
+ integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
+ double precision :: accu,hij
+ double precision :: norm
+ allocate (key_1(N_int,2),key_2(N_int,2))
+ do i = 1, n_act_orb
+  j = i           ! Diagonal part 
+  norm = 0.d0
+  accu = 0.d0
+  do k = 1,  n_det_ionic_on_given_atom(i)
+   norm += psi_coef_mono_ionic_on_given_atom(k,i) **2
+   do ii = 1, N_int
+    key_1(ii,1) = psi_det_mono_ionic_on_given_atom(ii,1,k,i)
+    key_1(ii,2) = psi_det_mono_ionic_on_given_atom(ii,2,k,i)
+   enddo 
+   do l = 1,  n_det_ionic_on_given_atom(j)
+    do jj = 1, N_int
+     key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j)
+     key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j)
+    enddo 
+    call i_H_j(key_1,key_2,N_int,hij)
+    accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
+   enddo
+  enddo
+  H_matrix_between_ionic_on_given_atom(i,j) = accu
+
+
+  do j = i+1, n_act_orb   ! Extra diagonal part 
+   accu = 0.d0
+   do k = 1,  n_det_ionic_on_given_atom(i)
+    do jj = 1, N_int
+     key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i)
+     key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i)
+    enddo 
+    do l = 1,  n_det_ionic_on_given_atom(j)
+     do jj = 1, N_int
+      key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j)
+      key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j)
+     enddo 
+     call i_H_j(key_1,key_2,N_int,hij)
+     accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
+    enddo
+   enddo
+   H_matrix_between_ionic_on_given_atom(i,j) = accu
+   H_matrix_between_ionic_on_given_atom(j,i) = accu
+  enddo
+ enddo
+ 
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom_and_others , (n_act_orb,min_number_ionic:max_number_ionic)]
+ implicit none
+ use bitmasks
+ BEGIN_DOC
+! Hamiltonian matrix elements between the various contracted functions 
+! that have a negative charge on a given active orbital 
+! and all the other fully contracted OVB structures 
+ END_DOC
+ integer :: i,j,k,l,jj,ii
+ integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
+ double precision :: accu,hij
+ double precision :: norm
+ allocate (key_1(N_int,2),key_2(N_int,2))
+ 
+ do i = 1, n_act_orb
+  do j = min_number_ionic,max_number_ionic
+   if(j==1)then
+    H_matrix_between_ionic_on_given_atom_and_others(i,j) = 0.d0
+   endif
+   accu = 0.d0
+   do k = 1,  n_det_ionic_on_given_atom(i)
+    do jj = 1, N_int
+     key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i)
+     key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i)
+    enddo 
+    do l = 1,  ionic_index(j,0)
+     do ii = 1, N_int
+      key_2(ii,1) = psi_det_ovb(ii,1,l,j)
+      key_2(ii,2) = psi_det_ovb(ii,2,l,j)
+     enddo 
+     call i_H_j(key_1,key_2,N_int,hij)
+     accu += psi_coef_ovb(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij
+    enddo
+   enddo
+   H_matrix_between_ionic_on_given_atom_and_others(i,j) = accu
+  enddo
+ enddo
+
+ print*,'H_matrix_between_ionic_on_given_atom_and_others'
+ print*,''
+ do i = 1, n_act_orb
+  write(*,'(I3,X,100(F16.7))'),H_matrix_between_ionic_on_given_atom_and_others(i,:)
+ enddo
+
+
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, n_det_ionic_on_given_atom, (n_act_orb)]
+&BEGIN_PROVIDER [double precision, normalization_factor_ionic_on_given_atom, (n_act_orb) ]
+&BEGIN_PROVIDER [double precision, psi_coef_mono_ionic_on_given_atom, (N_det_ref,n_act_orb) ]
+&BEGIN_PROVIDER [integer(bit_kind), psi_det_mono_ionic_on_given_atom, (N_int,2,N_det_ref,n_act_orb)]
+ implicit none
+ use bitmasks
+ BEGIN_DOC
+! number of determinants that are mono ionic with the negative charge 
+! on a given atom, normalization_factor, array of determinants,and coefficients
+ END_DOC
+ integer :: i,j,k,l
+ ionicity_level = 1
+ integer :: ionicity_level
+ logical :: doubly_occupied_array(n_act_orb)
+ n_det_ionic_on_given_atom = 0
+ normalization_factor_ionic_on_given_atom = 0.d0
+ do i = 1,  ionic_index(ionicity_level,0)
+  call give_index_of_doubly_occ_in_active_space(psi_det(1,1,ionic_index(ionicity_level,i)),doubly_occupied_array)
+  do j = 1, n_act_orb
+   if(doubly_occupied_array(j))then
+    n_det_ionic_on_given_atom(j) += 1
+    normalization_factor_ionic_on_given_atom(j) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2
+    do k = 1, N_int
+     psi_det_mono_ionic_on_given_atom(k,1,n_det_ionic_on_given_atom(j),j) = psi_det(k,1,ionic_index(ionicity_level,i))
+     psi_det_mono_ionic_on_given_atom(k,2,n_det_ionic_on_given_atom(j),j) = psi_det(k,2,ionic_index(ionicity_level,i))
+    enddo
+    psi_coef_mono_ionic_on_given_atom(n_det_ionic_on_given_atom(j),j) = psi_ref_coef_diagonalized(ionic_index(1,i),1) 
+   endif
+  enddo
+ enddo
+ integer :: i_count
+ i_count = 0
+ do j = 1, n_act_orb
+  i_count += n_det_ionic_on_given_atom(j)
+  normalization_factor_ionic_on_given_atom(j) = 1.d0/dsqrt(normalization_factor_ionic_on_given_atom(j))
+ enddo
+ if(i_count.ne.ionic_index(ionicity_level,0))then
+   print*,'PB with n_det_ionic_on_given_atom'
+   print*,'i_count = ',i_count
+   print*,'ionic_index(ionicity_level,0)',ionic_index(ionicity_level,0)
+   stop
+ endif
+ do j = 1, n_act_orb 
+  do i = 1, n_det_ionic_on_given_atom(j) 
+   psi_coef_mono_ionic_on_given_atom(i,j) =  psi_coef_mono_ionic_on_given_atom(i,j) * normalization_factor_ionic_on_given_atom(j)
+  enddo
+ enddo
+ 
+
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [integer(bit_kind), psi_det_ovb, (N_int,2,N_det_ref,min_number_ionic:max_number_ionic)]
+&BEGIN_PROVIDER [double precision, psi_coef_ovb, (N_det_ref,min_number_ionic:max_number_ionic)  ]
+ implicit none
+ BEGIN_DOC
+! Array of the determinants belonging to each ovb structures (neutral, mono ionic, bi ionic etc ...)
+! together with the arrays of coefficients
+ END_DOC
+ integer :: i,j,k,l
+ use bitmasks
+ integer :: ionicity_level,i_count
+ double precision :: accu
+ 
+ do ionicity_level = min_number_ionic,max_number_ionic
+  accu = 0.d0
+  do i = 1,  ionic_index(ionicity_level,0)
+   do j = 1, N_int
+    psi_det_ovb(j,1,i,ionicity_level)  = psi_det(j,1,ionic_index(ionicity_level,i))
+    psi_det_ovb(j,2,i,ionicity_level)  = psi_det(j,2,ionic_index(ionicity_level,i))
+   enddo
+   psi_coef_ovb(i,ionicity_level) = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) * normalization_factor_ionic(ionicity_level,1)
+   accu +=  psi_coef_ovb(i,ionicity_level)**2
+  enddo
+  accu = 1.d0/dsqrt(accu)
+  do i = 1,  ionic_index(ionicity_level,0)
+   psi_coef_ovb(i,ionicity_level) =  psi_coef_ovb(i,ionicity_level) * accu
+  enddo
+  accu = 0.d0
+  do i = 1,  ionic_index(ionicity_level,0)
+   accu += psi_coef_ovb(i,ionicity_level) **2
+  enddo
+ enddo
+ 
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, H_matrix_psi_det_ovb, (min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)]
+ implicit none
+ BEGIN_DOC
+! H matrix between the fully contracted OVB forms
+ END_DOC
+ integer :: i,j,k,l,jj,ii
+ integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:)
+ use bitmasks
+ double precision :: accu,hij
+ double precision :: norm
+ allocate (key_1(N_int,2),key_2(N_int,2))
+ do i = min_number_ionic,max_number_ionic
+  do j = min_number_ionic,max_number_ionic
+   accu = 0.d0
+   do k = 1, ionic_index(i,0)
+    do ii = 1, N_int
+     key_1(ii,1) = psi_det_ovb(ii,1,k,i)
+     key_1(ii,2) = psi_det_ovb(ii,2,k,i)
+    enddo
+    do l = 1, ionic_index(j,0)
+     do ii = 1, N_int
+      key_2(ii,1) = psi_det_ovb(ii,1,l,j)
+      key_2(ii,2) = psi_det_ovb(ii,2,l,j)
+     enddo
+     call i_H_j(key_1,key_2,N_int,hij)
+     accu += psi_coef_ovb(l,j) * psi_coef_ovb(k,i) * hij
+    enddo
+   enddo 
+   H_matrix_psi_det_ovb(i,j) = accu
+  enddo
+ enddo
+
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [integer, number_first_ionic_couples]
+&BEGIN_PROVIDER [logical , is_a_first_ionic_couple, (N_det_ref)]
+&BEGIN_PROVIDER [double precision, normalization_factor_special_first_ionic, (2)]
+ implicit none
+ BEGIN_DOC
+  ! Number of determinants belonging to the class of first ionic
+  ! AND that have a couple of positive/negative charge belonging 
+  ! to a couple of orbital couples_act_orb
+  ! If is_a_first_ionic_couple(i) = .True. then this determinant is a first ionic
+  ! and have a couple of positive/negative charge belonging
+  ! to a couple of orbital couples_act_orb
+  ! normalization factor (1) = 1/(sum c_i^2   .with. is_a_first_ionic_couple(i) = .True.)
+  ! normalization factor (2) = 1/(sum c_i^2   .with. is_a_first_ionic_couple(i) = .False.)
+ END_DOC
+ integer :: i,j
+ use bitmasks
+ number_first_ionic_couples = 0
+ integer :: ionicity_level
+ logical :: couples_out(0:n_couples_act_orb)
+ integer(bit_kind) :: key_tmp(N_int,2)
+ ionicity_level = 1
+ normalization_factor_special_first_ionic = 0.d0
+ do i = 1,  ionic_index(ionicity_level,0)
+  do j = 1, N_int
+   key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i))
+   key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i))
+  enddo
+  call doubly_occ_empty_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out)
+  if(couples_out(0))then
+   number_first_ionic_couples +=1
+   is_a_first_ionic_couple(i) = .True.
+   normalization_factor_special_first_ionic(1) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 
+  else 
+   is_a_first_ionic_couple(i) = .False.
+   normalization_factor_special_first_ionic(2) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2
+  endif
+ enddo
+ normalization_factor_special_first_ionic(1) = 1.d0/dsqrt(normalization_factor_special_first_ionic(1))
+ normalization_factor_special_first_ionic(2) = 1.d0/dsqrt(normalization_factor_special_first_ionic(2))
+ print*,'number_first_ionic_couples = ',number_first_ionic_couples
+ END_PROVIDER
+ 
+
+ BEGIN_PROVIDER [integer, number_neutral_no_hund_couples]
+&BEGIN_PROVIDER [logical , is_a_neutral_no_hund_couple, (N_det_ref)]
+&BEGIN_PROVIDER [double precision, normalization_factor_neutra_no_hund_couple, (2)]
+&BEGIN_PROVIDER [double precision, ratio_hund_no_hund ]
+ implicit none
+ BEGIN_DOC
+  ! Number of determinants belonging to the class of neutral determinants 
+  ! AND that have a couple of alpha beta electrons in couple of orbital couples_act_orb
+  ! If is_a_neutral_no_hund_couple(i) = .True. then this determinant is a neutral determinants 
+  ! and have a a couple of alpha beta electrons in couple of orbital couples_act_orb
+  ! normalization factor (1) = 1/sqrt(sum c_i^2   .with. is_a_neutral_no_hund_couple(i) = .True.)
+  ! normalization factor (2) = 1/sqrt(sum c_i^2   .with. is_a_neutral_no_hund_couple(i) = .False.)
+ END_DOC
+ integer :: i,j
+ use bitmasks
+ number_neutral_no_hund_couples = 0
+ integer :: ionicity_level
+ logical :: couples_out(0:n_couples_act_orb)
+ integer(bit_kind) :: key_tmp(N_int,2)
+ integer :: ifirst_hund,ifirst_no_hund
+ double precision :: coef_ref_hund,coef_ref_no_hund
+ ifirst_hund = 0
+ ifirst_no_hund = 0
+ ionicity_level = 0
+ normalization_factor_neutra_no_hund_couple = 0.d0
+ do i = 1,  ionic_index(ionicity_level,0)
+  do j = 1, N_int
+   key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i))
+   key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i))
+  enddo
+  call neutral_no_hund_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out)
+  if(couples_out(0))then
+   if(ifirst_no_hund == 0)then
+     coef_ref_no_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1)
+     ifirst_no_hund = 1
+   endif
+   number_neutral_no_hund_couples +=1
+   is_a_neutral_no_hund_couple(i) = .True.
+   normalization_factor_neutra_no_hund_couple(1) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2 
+  else 
+   if(ifirst_hund == 0)then
+     coef_ref_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1)
+     ifirst_hund = 1
+   endif
+   is_a_neutral_no_hund_couple(i) = .False.
+   normalization_factor_neutra_no_hund_couple(2) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2
+  endif
+ enddo
+ ratio_hund_no_hund = coef_ref_no_hund/coef_ref_hund
+
+ normalization_factor_neutra_no_hund_couple(1) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(1))
+ normalization_factor_neutra_no_hund_couple(2) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(2))
+ print*,'number_neutral_no_hund_couples = ',number_neutral_no_hund_couples
+ END_PROVIDER
+ 
+ BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic, (2,min_number_ionic:max_number_ionic,n_states)]
+&BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic_between_ionic, (2,2,n_states)]
+ BEGIN_DOC 
+ ! H_OVB_naked_first_ionic(1,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
+ ! and the contracted ith ionic form
+ ! if i == 1 not defined
+ ! H_OVB_naked_first_ionic(2,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False
+ ! and the contracted ith ionic form
+ ! if i == 1 not defined
+ ! H_OVB_naked_first_ionic_between_ionic(1,1) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
+ ! and the first ionic determinants belonging to is_a_first_ionic_couple = True
+ ! H_OVB_naked_first_ionic_between_ionic(1,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True
+ ! and the first ionic determinants belonging to is_a_first_ionic_couple = False
+ ! H_OVB_naked_first_ionic_between_ionic(2,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False
+ ! and the first ionic determinants belonging to is_a_first_ionic_couple = False
+ END_DOC
+ implicit none
+ integer :: i,j,istate,k,l
+ double precision :: accu_1,accu_2,hij
+  H_OVB_naked_first_ionic = 0.d0
+  H_OVB_naked_first_ionic_between_ionic = 0.d0
+  i = 1
+  do j = min_number_ionic,max_number_ionic
+   if(j==1)cycle
+   do istate = 1, N_states
+    accu_1 = 0.d0
+    accu_2 = 0.d0
+    do k = 1,  ionic_index(i,0)
+     if(is_a_first_ionic_couple(k))then   
+      do l = 1,  ionic_index(j,0)
+       hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
+       accu_1 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(1) * & 
+                 psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
+      enddo
+     else 
+      do l = 1,  ionic_index(j,0)
+       hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l))
+       accu_2 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(2) * & 
+                 psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
+      enddo
+     endif
+    enddo
+    H_OVB_naked_first_ionic(1,j,istate) = accu_1
+    H_OVB_naked_first_ionic(2,j,istate) = accu_2
+   enddo
+  enddo
+
+
+   do istate = 1, N_states
+    accu_1 = 0.d0
+    accu_2 = 0.d0
+    integer :: i_count
+    i_count = 0
+    do k = 1,  ionic_index(1,0)
+      do l = 1,  ionic_index(1,0)
+        hij = ref_hamiltonian_matrix(ionic_index(1,k),ionic_index(1,l))
+        accu_1 = hij * psi_ref_coef_diagonalized(ionic_index(1,k),istate) * psi_ref_coef_diagonalized(ionic_index(1,l),istate)
+        if(is_a_first_ionic_couple(k).and. is_a_first_ionic_couple(l))then   
+         H_OVB_naked_first_ionic_between_ionic(1,1,istate) += accu_1 *  normalization_factor_special_first_ionic(1) **2 
+        elseif(is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then
+         i_count += 1
+         H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 *  & 
+                                       normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2)
+!       elseif(is_a_first_ionic_couple(l).and. .not.is_a_first_ionic_couple(k))then
+!        i_count += 1
+!        H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 *  & 
+!                                      normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2)
+        elseif(.not.is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then
+         H_OVB_naked_first_ionic_between_ionic(2,2,istate) +=  accu_1 *  normalization_factor_special_first_ionic(2) **2
+        endif
+      enddo
+   enddo
+  enddo
+  print*,'i_count                       = ',i_count
+  print*,'number_first_ionic_couples**2 = ',ionic_index(1,0) * number_first_ionic_couples
+  
+ double precision :: convert_hartree_ev
+ convert_hartree_ev = 27.211399d0
+   print*,'Special H matrix'
+   do i = 1,2
+    write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic(i,:,1)
+   enddo
+
+   print*,'Special H matrix bis'
+   do i = 1,2
+    write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic_between_ionic(i,:,1)
+   enddo
+
+
+ END_PROVIDER 
+
diff --git a/plugins/OVB/print_ovb.irp.f b/plugins/OVB/print_ovb.irp.f
new file mode 100644
index 00000000..9e333c15
--- /dev/null
+++ b/plugins/OVB/print_ovb.irp.f
@@ -0,0 +1,27 @@
+program print_OVB
+ implicit none
+ read_wf = .True.
+ call provide_all
+
+end
+
+subroutine provide_all
+ implicit none
+ integer :: i,j,k,l,istate
+ do istate= 1, N_states
+ print*,'-------------------'
+ print*,'ISTATE = ',istate
+ print*,'-------------------'
+ print*,'CAS MATRIX         '
+ print*,''
+  do i = min_number_ionic,max_number_ionic
+   write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_naked(i,:,istate)
+  enddo
+  print*,''
+  print*,'-------------------'
+  print*,'-------------------'
+ enddo
+
+
+end
+
diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f
new file mode 100644
index 00000000..c1d88d2c
--- /dev/null
+++ b/plugins/Properties/hyperfine_constants.irp.f
@@ -0,0 +1,135 @@
+BEGIN_PROVIDER [double precision, spin_density_at_nucleous, (nucl_num)]
+ implicit none
+ BEGIN_DOC
+! value of the spin density at each nucleus 
+ END_DOC
+ integer :: i,j,k
+ do i = 1, nucl_num 
+  double precision :: r(3),accu,aos_array(ao_num)
+  accu = 0.d0
+  r(1:3) = nucl_coord(i,1:3)
+  call give_all_aos_at_r(r,aos_array)
+  do j = 1, ao_num
+   do k = 1, ao_num
+    accu += one_body_spin_density_ao(k,j) * aos_array(k) * aos_array(j)
+   enddo
+  enddo
+  spin_density_at_nucleous(i) = accu
+ enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, spin_density_at_nucleous_from_mo, (nucl_num)]
+&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_per_mo, (nucl_num,mo_tot_num)]
+ implicit none
+ BEGIN_DOC
+! value of the spin density at each nucleus 
+ END_DOC
+ integer :: i,j,k,l,m
+ do i = 1, nucl_num 
+  double precision :: r(3),accu,aos_array(ao_num)
+  double precision :: contrib
+  double precision :: mo_values(mo_tot_num)
+  accu = 0.d0
+  r(1:3) = nucl_coord(i,1:3)
+  call give_all_aos_at_r(r,aos_array)
+  spin_density_at_nucleous_from_mo(i) = 0.d0
+  do k = 1, mo_tot_num
+   mo_values(k) = 0.d0
+   do j = 1, ao_num
+    mo_values(k) += mo_coef(j,k) * aos_array(j)
+   enddo
+  enddo
+  do k = 1, mo_tot_num
+   spin_density_at_nucleous_contrib_per_mo(i,k) = 0.d0
+   do m = 1, mo_tot_num
+    contrib = one_body_spin_density_mo(k,m) * mo_values(k) * mo_values(m)
+    spin_density_at_nucleous_from_mo(i) += contrib
+    spin_density_at_nucleous_contrib_per_mo(i,k) += contrib
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo, (nucl_num,mo_tot_num,mo_tot_num)]
+&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo_test, (nucl_num)]
+ implicit none
+ BEGIN_DOC
+! value of the spin density at each nucleus 
+ END_DOC
+ integer :: i,j,k,l,m
+ spin_density_at_nucleous_contrib_mo_test = 0.d0
+ do i = 1, nucl_num 
+  double precision :: r(3),accu,aos_array(ao_num)
+  double precision :: c_i1,c_j1  
+  r(1:3) = nucl_coord(i,1:3)
+  call give_all_aos_at_r(r,aos_array)
+  do k = 1, mo_tot_num
+   do m = 1, mo_tot_num
+    accu = 0.d0
+       do j = 1, ao_num
+        c_i1 = mo_coef(j,k)
+        do l = 1, ao_num
+         c_j1 = c_i1*mo_coef(l,m)
+         accu += one_body_spin_density_mo(k,m) * aos_array(l) * aos_array(j) * c_j1
+        enddo
+       enddo
+    spin_density_at_nucleous_contrib_mo(i,k,m) = accu
+    spin_density_at_nucleous_contrib_mo_test(i) += accu
+   enddo 
+  enddo
+ enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, conversion_factor_mhz_hcc, (100)]
+&BEGIN_PROVIDER [double precision, conversion_factor_gauss_hcc, (100)]
+&BEGIN_PROVIDER [double precision, conversion_factor_cm_1_hcc, (100)]
+ BEGIN_DOC
+! Conversion factor for the calculation of the hcc, according to the nuclear charge 
+ END_DOC
+
+ conversion_factor_mhz_hcc =0.d0 
+ conversion_factor_mhz_hcc =0.d0 
+ conversion_factor_mhz_hcc =0.d0 
+
+
+ ! hydrogen 
+ conversion_factor_mhz_hcc(1) = 4469.84692227102460d0
+ conversion_factor_gauss_hcc(1) = 1594.95296390862904d0
+ conversion_factor_cm_1_hcc(1) = 1490.98044430157870d0
+
+ ! Li       
+ conversion_factor_mhz_hcc(3) = 1737.2746512855997d0
+ conversion_factor_gauss_hcc(3) = 619.9027742370165d0
+ conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
+
+ ! carbon
+ conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
+ conversion_factor_gauss_hcc(6) = 401.136570647523058d0
+ conversion_factor_cm_1_hcc(6) = 374.987097339830086d0
+
+ ! nitrogen
+ conversion_factor_mhz_hcc(7) = 323.102093833793390d0
+ conversion_factor_gauss_hcc(7) = 115.290892768082614d0
+ conversion_factor_cm_1_hcc(7) = 107.775257586297698d0
+
+ ! Oxygen
+ conversion_factor_mhz_hcc(8) = -606.1958551736545d0
+ conversion_factor_gauss_hcc(8) = -216.30574771560407d0
+ conversion_factor_cm_1_hcc(8) = -202.20517197179822d0
+ 
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, iso_hcc_mhz, (nucl_num)]
+&BEGIN_PROVIDER [double precision, iso_hcc_gauss, (nucl_num)]
+&BEGIN_PROVIDER [double precision, iso_hcc_cm_1, (nucl_num)]
+ BEGIN_DOC
+! isotropic hyperfine coupling constants among the various atoms
+ END_DOC
+ integer :: i
+ do i = 1, nucl_num
+   iso_hcc_mhz(i)   = conversion_factor_mhz_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)  !* 0.5d0
+   iso_hcc_gauss(i) = conversion_factor_gauss_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)!* 0.5d0
+   iso_hcc_cm_1(i)  = conversion_factor_cm_1_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !*0.5d0
+ enddo
+
+END_PROVIDER 
diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f
new file mode 100644
index 00000000..d56c9a44
--- /dev/null
+++ b/plugins/Properties/mulliken.irp.f
@@ -0,0 +1,107 @@
+
+BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)]
+ implicit none
+ integer :: i,j
+ BEGIN_DOC
+! spin population on the ao basis :
+! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
+ END_DOC
+ spin_population = 0.d0
+ do i = 1, ao_num
+  do j = 1, ao_num
+   spin_population(j,i) = one_body_spin_density_ao(i,j) * ao_overlap(i,j)
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)]
+ implicit none
+ integer :: i
+ double precision :: accu
+ spin_population_angular_momentum = 0.d0
+ do i = 1,  ao_num
+  spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i)
+ enddo
+
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [double precision, spin_gross_orbital_product, (ao_num)]
+ implicit none
+ spin_gross_orbital_product = 0.d0
+ integer :: i,j
+ BEGIN_DOC
+! gross orbital product for the spin population
+ END_DOC
+ do i = 1, ao_num
+  do j = 1, ao_num
+   spin_gross_orbital_product(i) += spin_population(j,i)
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, mulliken_spin_densities, (nucl_num)]
+ implicit none
+ integer :: i,j
+ BEGIN_DOC
+!ATOMIC SPIN POPULATION (ALPHA MINUS BETA)
+ END_DOC
+ mulliken_spin_densities = 0.d0
+ do i = 1, ao_num
+  mulliken_spin_densities(ao_nucl(i)) += spin_gross_orbital_product(i)
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, electronic_population_alpha, (ao_num_align,ao_num)]
+&BEGIN_PROVIDER [double precision, electronic_population_beta, (ao_num_align,ao_num)]
+ implicit none
+ integer :: i,j
+ BEGIN_DOC
+! spin population on the ao basis :
+! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
+ END_DOC
+ electronic_population_alpha = 0.d0
+ electronic_population_beta = 0.d0
+ do i = 1, ao_num
+  do j = 1, ao_num
+   electronic_population_alpha(j,i) = one_body_dm_ao_alpha(i,j) * ao_overlap(i,j)
+   electronic_population_beta(j,i) = one_body_dm_ao_beta(i,j) * ao_overlap(i,j)
+  enddo
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, gross_orbital_product_alpha, (ao_num)]
+&BEGIN_PROVIDER [double precision, gross_orbital_product_beta, (ao_num)]
+ implicit none
+ spin_gross_orbital_product = 0.d0
+ integer :: i,j
+ BEGIN_DOC
+! gross orbital product 
+ END_DOC
+ do i = 1, ao_num
+  do j = 1, ao_num
+   gross_orbital_product_alpha(i) += electronic_population_alpha(j,i)
+   gross_orbital_product_beta(i) += electronic_population_beta(j,i)
+  enddo
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, mulliken_densities_alpha, (nucl_num)]
+&BEGIN_PROVIDER [double precision, mulliken_densities_beta, (nucl_num)]
+ implicit none
+ integer :: i,j
+ BEGIN_DOC
+!
+ END_DOC
+ mulliken_densities_alpha = 0.d0
+ mulliken_densities_beta = 0.d0
+ do i = 1, ao_num
+  mulliken_densities_alpha(ao_nucl(i)) += gross_orbital_product_alpha(i)
+  mulliken_densities_beta(ao_nucl(i)) += gross_orbital_product_beta(i)
+ enddo
+
+END_PROVIDER
diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f
new file mode 100644
index 00000000..f0091e1e
--- /dev/null
+++ b/plugins/Properties/print_hcc.irp.f
@@ -0,0 +1,17 @@
+program print_hcc
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call test
+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
+
diff --git a/plugins/Properties/print_mulliken.irp.f b/plugins/Properties/print_mulliken.irp.f
new file mode 100644
index 00000000..100c8556
--- /dev/null
+++ b/plugins/Properties/print_mulliken.irp.f
@@ -0,0 +1,35 @@
+program print_mulliken
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ print*,'Mulliken spin densities'
+
+ 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
+
diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f
index f9cf1303..fb45b13d 100644
--- a/plugins/Psiref_Utils/psi_ref_utils.irp.f
+++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f
@@ -125,7 +125,7 @@ BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)]
   enddo
 END_PROVIDER
 
- BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)]
+ BEGIN_PROVIDER [double precision, psi_ref_coef_diagonalized, (N_det_ref,N_states)]
 &BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)]
  implicit none
  integer :: i,j
@@ -137,9 +137,11 @@ END_PROVIDER
   do i = 1, N_states
    psi_ref_energy_diagonalized(i) = eigenvalues(i)
    do j = 1, N_det_ref
-    psi_coef_ref_diagonalized(j,i) = eigenvectors(j,i)
+    psi_ref_coef_diagonalized(j,i) = eigenvectors(j,i)
    enddo
   enddo
+  deallocate (eigenvectors)
+  deallocate (eigenvalues)
 
 
  END_PROVIDER
@@ -264,3 +266,18 @@ integer function get_index_in_psi_ref_sorted_bit(key,Nint)
 
 end
 
+BEGIN_PROVIDER [double precision, ref_hamiltonian_matrix, (n_det_ref,n_det_ref)]
+ BEGIN_DOC
+ ! H matrix in the Reference space
+ END_DOC
+ implicit none
+ integer :: i,j
+ double precision :: hij
+ do i = 1, N_det_ref
+  do j = 1, N_det_ref
+   call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
+   ref_hamiltonian_matrix(i,j) = hij
+  enddo
+ enddo
+END_PROVIDER
+
diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py
index b52e4445..f83d8b41 100755
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@@ -156,11 +156,11 @@ class H_apply(object):
   def filter_only_1h1p(self):
     self["filter_only_1h1p_single"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h1p(hole).eq..False.) cycle
+     if (is_a_1h1p(hole).eqv..False.) cycle
     """
     self["filter_only_1h1p_double"] = """
 !    ! DIR$ FORCEINLINE
-     if (is_a_1h1p(key).eq..False.) cycle
+     if (is_a_1h1p(key).eqv..False.) cycle
     """
 
 
diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f
index 341d1453..8c2db90e 100644
--- a/src/AO_Basis/aos.irp.f
+++ b/src/AO_Basis/aos.irp.f
@@ -146,3 +146,30 @@ integer function ao_power_index(nx,ny,nz)
   ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
 end
 
+ BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
+&BEGIN_PROVIDER [ integer, ao_l_max ]
+&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ]
+ implicit none
+ BEGIN_DOC
+! ao_l = l value of the AO: a+b+c in x^a y^b z^c
+ END_DOC
+ integer :: i
+ do i=1,ao_num
+   ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
+   ao_l_char(i) = l_to_charater(ao_l(i))
+ enddo
+ ao_l_max = maxval(ao_l)
+END_PROVIDER
+
+BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
+ BEGIN_DOC
+ ! character corresponding to the "L" value of an AO orbital
+ END_DOC
+ implicit none
+ l_to_charater(0)='S'
+ l_to_charater(1)='P'
+ l_to_charater(2)='D'
+ l_to_charater(3)='F'
+ l_to_charater(4)='G'
+END_PROVIDER
+
diff --git a/src/AO_Basis/aos_value.irp.f b/src/AO_Basis/aos_value.irp.f
new file mode 100644
index 00000000..a531ce50
--- /dev/null
+++ b/src/AO_Basis/aos_value.irp.f
@@ -0,0 +1,48 @@
+double precision function ao_value(i,r)
+ implicit none
+ BEGIN_DOC
+! return the value of the ith ao at point r
+ END_DOC
+ double precision, intent(in) :: r(3)
+ integer, intent(in) :: i
+
+ integer :: m,num_ao
+ double precision :: center_ao(3)
+ double precision :: beta
+ integer :: power_ao(3) 
+ num_ao = ao_nucl(i)
+ power_ao(1:3)= ao_power(i,1:3)
+ center_ao(1:3) = nucl_coord(num_ao,1:3)
+ double precision :: accu,dx,dy,dz,r2
+ dx = (r(1) - center_ao(1)) 
+ dy = (r(2) - center_ao(2)) 
+ dz = (r(3) - center_ao(3)) 
+ r2 = dx*dx + dy*dy + dz*dz
+ dx = dx**power_ao(1)
+ dy = dy**power_ao(2)
+ dz = dz**power_ao(3)
+
+ accu = 0.d0
+ do m=1,ao_prim_num(i)
+   beta = ao_expo_ordered_transp(m,i)
+   accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) 
+ enddo
+ ao_value = accu * dx * dy * dz
+
+end
+
+subroutine give_all_aos_at_r(r,aos_array)
+ implicit none
+ BEGIN_dOC
+! gives the values of aos at a given point r
+ END_DOC
+ double precision, intent(in) :: r(3)
+ double precision, intent(out) :: aos_array(ao_num)
+ integer :: i
+ double precision :: ao_value
+ do i = 1, ao_num 
+  aos_array(i) = ao_value(i,r)
+ enddo
+
+
+end
diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f
index 70a84a42..50a95312 100644
--- a/src/Bitmask/bitmasks.irp.f
+++ b/src/Bitmask/bitmasks.irp.f
@@ -289,7 +289,12 @@ END_PROVIDER
 &BEGIN_PROVIDER [ integer, n_virt_orb ]
  implicit none
  BEGIN_DOC
- ! Bitmasks for the inactive orbitals that are excited in post CAS method
+ ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited 
+ ! in post CAS methods
+ ! n_inact_orb   : Number of inactive orbitals
+ ! virt_bitmask  : Bitmaks of vritual orbitals which are supposed to be recieve electrons 
+ ! in post CAS methods
+ ! n_virt_orb    : Number of virtual orbitals
  END_DOC
  logical                        :: exists
  integer                        :: j,i
@@ -327,8 +332,14 @@ END_PROVIDER
 
 
 
- BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)]
+  BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)]
  &BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)]
+ BEGIN_DOC
+ ! list_inact : List of the inactive orbitals which are supposed to be doubly excited 
+ ! in post CAS methods
+ ! list_virt  : List of vritual orbitals which are supposed to be recieve electrons 
+ ! in post CAS methods
+ END_DOC
  implicit none
  integer :: occ_inact(N_int*bit_kind_size)
  integer :: itest,i
@@ -348,6 +359,21 @@ END_PROVIDER
 
  END_PROVIDER 
 
+ BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)]
+ implicit none
+ BEGIN_DOC
+ ! Reunion of the inactive, active and virtual bitmasks
+ END_DOC
+ integer :: i,j
+ do i = 1, N_int
+  reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1))
+  reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2))
+ enddo
+ END_PROVIDER
+
+
+
+
  BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)]
  implicit none
  BEGIN_DOC
@@ -376,7 +402,7 @@ END_PROVIDER
  BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)]
  implicit none
  BEGIN_DOC
- ! Reunion of the inactive, active and virtual bitmasks
+ ! Bitmask of the core orbitals that are never excited in post CAS method
  END_DOC
  integer :: i,j
  do i = 1, N_int
@@ -385,6 +411,35 @@ END_PROVIDER
  enddo
  END_PROVIDER 
 
+ BEGIN_PROVIDER [integer, list_core, (n_core_orb)]
+ BEGIN_DOC
+ ! List of the core orbitals that are never excited in post CAS method
+ END_DOC
+ implicit none
+ integer :: occ_core(N_int*bit_kind_size)
+ integer :: itest,i
+ occ_core = 0
+ call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int)
+ ASSERT(itest==n_core_orb)
+ do i = 1, n_core_orb
+  list_core(i) = occ_core(i)
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, n_core_orb ]
+ implicit none
+ BEGIN_DOC
+ ! Number of core orbitals that are never excited in post CAS method
+ END_DOC
+ logical                        :: exists
+ integer                        :: j,i
+ integer :: i_hole,i_part,i_gen
+
+ n_core_orb = 0
+ do j = 1, N_int
+  n_core_orb += popcnt(core_bitmask(j,1))
+ enddo
+ END_PROVIDER 
 
 
 BEGIN_PROVIDER [ integer, i_bitmask_gen ]
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index e5d243f4..62d09381 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -206,3 +206,54 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
  state_average_weight = 1.d0/dble(N_states)
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, one_body_spin_density_ao, (ao_num_align,ao_num) ]
+ BEGIN_DOC
+! one body spin density matrix on the AO basis : rho_AO(alpha) - rho_AO(beta)
+ END_DOC
+ implicit none
+ integer :: i,j,k,l
+ double precision :: dm_mo
+
+ one_body_spin_density_ao = 0.d0
+ do k = 1, ao_num
+  do l = 1, ao_num
+   do i = 1, mo_tot_num
+    do j = 1, mo_tot_num
+     dm_mo = one_body_spin_density_mo(j,i)
+!    if(dabs(dm_mo).le.1.d-10)cycle
+     one_body_spin_density_ao(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
+
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ]
+ BEGIN_DOC
+! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta)
+ END_DOC
+ implicit none
+ integer :: i,j,k,l
+ double precision :: dm_mo
+
+ one_body_spin_density_ao = 0.d0
+ do k = 1, ao_num
+  do l = 1, ao_num
+   do i = 1, mo_tot_num
+    do j = 1, mo_tot_num
+     dm_mo = one_body_dm_mo_alpha(j,i)
+!    if(dabs(dm_mo).le.1.d-10)cycle
+     one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
+     one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
+
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
diff --git a/src/Determinants/usefull_for_ovb.irp.f b/src/Determinants/usefull_for_ovb.irp.f
new file mode 100644
index 00000000..7b89897b
--- /dev/null
+++ b/src/Determinants/usefull_for_ovb.irp.f
@@ -0,0 +1,283 @@
+
+integer function n_open_shell(det_in,nint)
+ implicit none
+ use bitmasks
+ integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer :: i
+ n_open_shell = 0
+ do i=1,Nint
+   n_open_shell += popcnt(iand(xor(det_in(i,1),det_in(i,2)),det_in(i,1)))
+ enddo
+end
+
+integer function n_closed_shell(det_in,nint)
+ implicit none
+ use bitmasks
+ integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer :: i
+ n_closed_shell = 0
+ do i=1,Nint
+   n_closed_shell += popcnt(iand(det_in(i,1),det_in(i,2)))
+ enddo
+end
+
+integer function n_closed_shell_cas(det_in,nint)
+ implicit none
+ use bitmasks
+ integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer(bit_kind) :: det_tmp(nint,2)
+ integer :: i
+ n_closed_shell_cas = 0
+ do i=1,Nint
+   det_tmp(i,1) = xor(det_in(i,1),reunion_of_core_inact_bitmask(i,1))
+   det_tmp(i,2) = xor(det_in(i,2),reunion_of_core_inact_bitmask(i,2))
+ enddo
+!call debug_det(det_tmp,nint)
+ do i=1,Nint
+   n_closed_shell_cas += popcnt(iand(det_tmp(i,1),det_tmp(i,2)))
+ enddo
+end
+
+subroutine doubly_occ_empty_in_couple(det_in,n_couples,couples,couples_out)
+ implicit none
+  use bitmasks
+ integer, intent(in) :: n_couples,couples(n_couples,2)
+ integer(bit_kind),intent(in) :: det_in(N_int,2)
+ logical, intent(out) :: couples_out(0:n_couples)
+ integer(bit_kind) :: det_tmp(N_int)
+ integer(bit_kind) :: det_tmp_bis(N_int)
+ BEGIN_DOC
+ ! n_couples is the number of couples of orbitals to be checked
+ ! couples(i,1) = first orbital of the ith couple
+ ! couples(i,2) = second orbital of the ith couple
+ ! returns the array couples_out
+ ! couples_out(i) = .True. if det_in contains 
+ !                   an orbital empty in the ith couple  AND
+ !                   an orbital doubly occupied in the ith couple
+ END_DOC
+ integer :: i,j,k,l
+ 
+ ! det_tmp tells you if the orbitals are occupied or not
+ do j = 1, N_int
+  det_tmp(j) = ior(det_in(j,1),det_in(j,2))
+ enddo
+ 
+ couples_out(0) = .False.
+ do i = 1, n_couples
+  do j = 1, N_int
+   det_tmp_bis(j) = 0_bit_kind
+  enddo
+  call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
+  call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
+  ! det_tmp is zero except for the two orbitals of the couple
+  integer :: i_count
+  i_count = 0
+  do j = 1, N_int
+    i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j)))   ! check if the two orbitals are both occupied 
+  enddo
+  if(i_count .ne. 1)then
+   couples_out(i) = .False.  
+   cycle  
+  endif
+
+  ! test if orbital there are two electrons or not
+  i_count = 0
+  do j = 1, N_int
+   i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
+  enddo
+  if(i_count.ne.1)then
+   couples_out(i) = .False.
+  else  
+   couples_out(i) = .True.
+   couples_out(0) = .True.
+  endif
+ enddo
+end
+
+subroutine give_index_of_doubly_occ_in_active_space(det_in,doubly_occupied_array)
+ implicit none
+  use bitmasks
+ integer(bit_kind), intent(in) :: det_in(N_int,2)
+ logical, intent(out) :: doubly_occupied_array(n_act_orb)
+ integer(bit_kind) :: det_tmp(N_int)
+ integer(bit_kind) :: det_tmp_bis(N_int)
+ BEGIN_DOC
+ END_DOC
+ integer :: i,j,k,l
+ 
+ ! det_tmp tells you if the orbitals are occupied or not
+ do j = 1, N_int
+  det_tmp(j) = ior(det_in(j,1),det_in(j,2))
+ enddo
+ 
+ do i = 1, n_act_orb 
+  do j = 1, N_int
+   det_tmp_bis(j) = 0_bit_kind
+  enddo
+  i_bite = list_act(i)
+  call set_bit_to_integer(i_bite,det_tmp_bis,N_int) ! act orb
+  ! det_tmp is zero except for the orbital "ith" active orbital
+  integer :: i_count,i_bite
+
+  ! test if orbital there are two electrons or not
+  i_count = 0
+  do j = 1, N_int
+   i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
+  enddo
+  if(i_count.ne.1)then
+   doubly_occupied_array(i) = .False.
+  else  
+   doubly_occupied_array(i) = .True.
+  endif
+ enddo
+end
+
+subroutine doubly_occ_empty_in_couple_and_no_hund_elsewhere(det_in,n_couple_no_hund,couple_ion,couple_no_hund,is_ok)
+ implicit none
+  use bitmasks
+ integer, intent(in) :: n_couple_no_hund,couple_ion(2),couple_no_hund(n_couple_no_hund,2)
+ integer(bit_kind),intent(in) :: det_in(N_int,2)
+ logical, intent(out) :: is_ok
+ integer(bit_kind) :: det_tmp(N_int)
+ integer(bit_kind) :: det_tmp_bis(N_int)
+ BEGIN_DOC
+ ! n_couples is the number of couples of orbitals to be checked
+ ! couples(i,1) = first orbital of the ith couple
+ ! couples(i,2) = second orbital of the ith couple
+ ! returns the array couples_out
+ ! couples_out(i) = .True. if det_in contains 
+ !                   an orbital empty in the ith couple  AND
+ !                   an orbital doubly occupied in the ith couple
+ END_DOC
+ integer :: i,j,k,l
+ 
+ ! det_tmp tells you if the orbitals are occupied or not
+ do j = 1, N_int
+  det_tmp(j) = ior(det_in(j,1),det_in(j,2))
+ enddo
+ 
+  is_ok = .False.
+  do j = 1, N_int
+   det_tmp_bis(j) = 0_bit_kind
+  enddo
+  call set_bit_to_integer(couple_ion(1),det_tmp_bis,N_int) ! first orb
+  call set_bit_to_integer(couple_ion(2),det_tmp_bis,N_int) ! second orb
+  ! det_tmp is zero except for the two orbitals of the couple
+  integer :: i_count
+  i_count = 0
+  do j = 1, N_int
+    i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j)))   ! check if the two orbitals are both occupied 
+  enddo
+  if(i_count .ne. 1)then
+   is_ok = .False.  
+   return
+  endif
+
+  ! test if orbital there are two electrons or not
+  i_count = 0
+  do j = 1, N_int
+   i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
+  enddo
+  if(i_count.ne.1)then
+   is_ok = .False.  
+   return
+  else  
+   do i = 1, n_couple_no_hund
+    do j = 1, N_int
+     det_tmp_bis(j) = 0_bit_kind
+    enddo
+    call set_bit_to_integer(couple_no_hund (i,1),det_tmp_bis,N_int) ! first orb
+    call set_bit_to_integer(couple_no_hund (i,2),det_tmp_bis,N_int) ! second orb
+    ! det_tmp_bis is zero except for the two orbitals of the couple
+    i_count = 0
+    do j = 1, N_int
+      i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j)))   ! check if the two orbitals are both occupied 
+    enddo
+    if(i_count .ne. 2)then
+     is_ok = .False.  
+     return
+    endif
+    ! test if orbital there are one alpha and one beta
+    integer :: i_count_alpha,i_count_beta
+    i_count_alpha = 0
+    i_count_beta = 0
+    do j = 1, N_int
+     i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
+     i_count_beta  += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
+    enddo
+    if(i_count_alpha==1.and.i_count_beta==1)then
+     is_ok = .True.
+    else
+     is_ok = .False.
+     return
+    endif
+   enddo
+   is_ok = .True.
+  endif
+end
+
+
+subroutine neutral_no_hund_in_couple(det_in,n_couples,couples,couples_out)
+ implicit none
+  use bitmasks
+ integer, intent(in) :: n_couples,couples(n_couples,2)
+ integer(bit_kind),intent(in) :: det_in(N_int,2)
+ logical, intent(out) :: couples_out(0:n_couples)
+ integer(bit_kind) :: det_tmp(N_int)
+ integer(bit_kind) :: det_tmp_bis(N_int)
+ BEGIN_DOC
+ ! n_couples is the number of couples of orbitals to be checked
+ ! couples(i,1) = first orbital of the ith couple
+ ! couples(i,2) = second orbital of the ith couple
+ ! returns the array couples_out
+ ! couples_out(i) = .True. if det_in contains 
+ !                   an orbital empty in the ith couple  AND
+ !                   an orbital doubly occupied in the ith couple
+ END_DOC
+ integer :: i,j,k,l
+ 
+ ! det_tmp tells you if the orbitals are occupied or not
+ do j = 1, N_int
+  det_tmp(j) = ior(det_in(j,1),det_in(j,2))
+ enddo
+ 
+ couples_out(0) = .True.
+ do i = 1, n_couples
+  do j = 1, N_int
+   det_tmp_bis(j) = 0_bit_kind
+  enddo
+  call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
+  call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
+  ! det_tmp_bis is zero except for the two orbitals of the couple
+  integer :: i_count
+  i_count = 0
+  do j = 1, N_int
+    i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j)))   ! check if the two orbitals are both occupied 
+  enddo
+  if(i_count .ne. 2)then
+   couples_out(i) = .False.  
+   cycle  
+  endif
+
+  ! test if orbital there are one alpha and one beta
+  integer :: i_count_alpha,i_count_beta
+  i_count_alpha = 0
+  i_count_beta = 0
+  do j = 1, N_int
+   i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
+   i_count_beta  += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
+  enddo
+  if(i_count_alpha==1.and.i_count_beta==1)then
+   couples_out(i) = .True.
+  else
+   couples_out(i) = .False.
+  endif
+ enddo
+ do i = 1, n_couples
+  if(.not.couples_out(i))then
+   couples_out(0) = .False.
+  endif
+ enddo
+end
+
+
diff --git a/src/Integrals_Bielec/.gitignore b/src/Integrals_Bielec/.gitignore
index ad6b465d..1d52a821 100644
--- a/src/Integrals_Bielec/.gitignore
+++ b/src/Integrals_Bielec/.gitignore
@@ -17,5 +17,4 @@ ZMQ
 ezfio_interface.irp.f
 irpf90.make
 irpf90_entities
-tags
-test_integrals
\ No newline at end of file
+tags
\ No newline at end of file
diff --git a/src/Integrals_Monoelec/.gitignore b/src/Integrals_Monoelec/.gitignore
index 577068de..e8bd9b05 100644
--- a/src/Integrals_Monoelec/.gitignore
+++ b/src/Integrals_Monoelec/.gitignore
@@ -12,9 +12,7 @@ Makefile.depend
 Nuclei
 Pseudo
 Utils
-check_orthonormality
 ezfio_interface.irp.f
 irpf90.make
 irpf90_entities
-save_ortho_mos
 tags
\ No newline at end of file
diff --git a/src/MOGuess/.gitignore b/src/MOGuess/.gitignore
index e9ba5cf5..797574f4 100644
--- a/src/MOGuess/.gitignore
+++ b/src/MOGuess/.gitignore
@@ -4,7 +4,6 @@
 AO_Basis
 Electrons
 Ezfio_files
-H_CORE_guess
 IRPF90_man
 IRPF90_temp
 Integrals_Monoelec
@@ -15,8 +14,6 @@ Nuclei
 Pseudo
 Utils
 ezfio_interface.irp.f
-guess_overlap
 irpf90.make
 irpf90_entities
-tags
-truncate_mos
\ No newline at end of file
+tags
\ No newline at end of file