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/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_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/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/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/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