diff --git a/plugins/local/ao_extra_basis/EZFIO.cfg b/plugins/local/ao_extra_basis/EZFIO.cfg
index 8b3a8667..8841c194 100644
--- a/plugins/local/ao_extra_basis/EZFIO.cfg
+++ b/plugins/local/ao_extra_basis/EZFIO.cfg
@@ -3,6 +3,12 @@ type: character*(256)
 doc: Name of the |ao_extra| basis set
 interface: ezfio
 
+[ao_extra_only_1s]
+type: logical
+doc: If |true|, you know that the additional AO basis is built only with 1s functions 
+interface: ezfio, provider
+default: true
+
 [ao_extra_num]
 type: integer
 doc: Number of |ao_extras|
@@ -85,3 +91,14 @@ doc: phase shift for each primitive GTOs |ao_extra|
 size: (3,ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
 interface: ezfio, provider
 
+[ao_extra_one_e_dm]
+type: double precision
+doc: reduced density matrix on the ao extra basis 
+size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_num)
+interface: ezfio, provider
+
+[ao_extra_center]
+type: double precision
+doc: shift with which the atoms are shifted to mimick p functions
+interface: ezfio
+
diff --git a/plugins/local/ao_extra_basis/LiH.xyz b/plugins/local/ao_extra_basis/LiH.xyz
new file mode 100644
index 00000000..03431803
--- /dev/null
+++ b/plugins/local/ao_extra_basis/LiH.xyz
@@ -0,0 +1,4 @@
+2
+
+H 0. 0. 0.
+Li 0. 0. 1.0
diff --git a/plugins/local/ao_extra_basis/README.rst b/plugins/local/ao_extra_basis/README.rst
index 3d1dc813..f60d71c0 100644
--- a/plugins/local/ao_extra_basis/README.rst
+++ b/plugins/local/ao_extra_basis/README.rst
@@ -2,3 +2,14 @@
 extra_basis
 ===========
 
+Plugin to handle an extra basis, which is attached to the extra_nuclei.
+It is essentially a duplication of all important quantities (coefficients, exponents and so on) of the usual |AO| basis. 
+
+An interesting feature is the possibility to fit any basis made at most with "p" functions onto a purely "s" basis. 
+This is done with the various scripts here: 
+ 
+ - qp_fit_1s_basis : script that creates an |EZFIO| folder corresponding to an .xyz file and a basis fitted with only "s" functions
+ - qp_add_extra_fit_system : script that takes as input an |EZFIO| folder and an .xyz file and add an extra_basis and extra_nuclei with a 1s fitted basis
+ 
+Ex: 
+qp_add_extra_fit_system LiH.ezfio/ h2o.xyz # takes the EZFIO folder "LiH.ezfio" and creates all necessary additional basis and nuclei based on h2o.xyz, but only with 1s functions. 
diff --git a/plugins/local/ao_extra_basis/fit_1s_basis.irp.f b/plugins/local/ao_extra_basis/fit_1s_basis.irp.f
index 3be2dbff..ef09d5b2 100644
--- a/plugins/local/ao_extra_basis/fit_1s_basis.irp.f
+++ b/plugins/local/ao_extra_basis/fit_1s_basis.irp.f
@@ -10,11 +10,28 @@ program fit_1s_basis
  print*,'Writting the results in the extra_nuclei and ao_extra_basis folders of EZFIO'
  print*,'New number of atomic functions : '
  print*,'n_func_tot = ',n_func_tot
+
+ print*,'extra_fictious_nucl = ',extra_fictious_nucl
+ do i = 1, extra_fictious_nucl
+  print*,list_fict_nucl(i)
+ enddo
+ print*,''
+ print*,''
+ do i = 1, nucl_num
+  print*,list_real_nucl(i)
+ enddo
  call ezfio_set_extra_nuclei_extra_nucl_num(new_nucl_num)
+ call ezfio_set_extra_nuclei_extra_nucl_fictious_num(extra_fictious_nucl)
+ call ezfio_set_extra_nuclei_extra_nucl_real_num(nucl_num)
+ call ezfio_set_extra_nuclei_extra_nucl_fictious_list(list_fict_nucl)
+ call ezfio_set_extra_nuclei_extra_nucl_real_list(list_real_nucl)
+ call ezfio_set_extra_nuclei_extra_nucl_real_fictious_list(extra_nucl_real_fictious_list_prov)
  call ezfio_set_extra_nuclei_extra_nucl_charge(new_nucl_charge_1s)
  call ezfio_set_extra_nuclei_extra_nucl_coord(new_nucl_coord_1s)
  call ezfio_set_extra_nuclei_extra_nucl_label(new_nucl_label_1s)
+!
  call ezfio_set_ao_extra_basis_ao_extra_num(n_func_tot)
+ call ezfio_set_ao_extra_basis_ao_extra_center(ao_extra_center)
  call ezfio_set_ao_extra_basis_ao_extra_nucl(new_ao_nucl_1s)
  call ezfio_set_ao_extra_basis_ao_extra_prim_num(new_ao_prim_num_1s)
  call ezfio_set_ao_extra_basis_ao_extra_coef(new_ao_coef_1s)
diff --git a/plugins/local/ao_extra_basis/h2o.xyz b/plugins/local/ao_extra_basis/h2o.xyz
new file mode 100644
index 00000000..d3928214
--- /dev/null
+++ b/plugins/local/ao_extra_basis/h2o.xyz
@@ -0,0 +1,7 @@
+3
+
+O      0.000000    -0.399441  3.000000   
+H      0.761232     0.199721  3.000000   
+H     -0.761232     0.199721  3.000000   
+
+
diff --git a/plugins/local/ao_extra_basis/install b/plugins/local/ao_extra_basis/install
index ca4a1643..e7e668e1 100755
--- a/plugins/local/ao_extra_basis/install
+++ b/plugins/local/ao_extra_basis/install
@@ -12,7 +12,7 @@ fi
 currdir=${PWD}
 
 # list of the scripts to be used by the module 
-scripts_list="qp_copy_extra_basis"
+scripts_list="qp_copy_extra_basis qp_add_extra_fit_system qp_copy_extra_basis_to_usual_basis qp_fit_1s_basis"
 
 # Make a symbolic link for all scripts to be used in the ${QP_ROOT}/scripts/
 # directory. 
diff --git a/plugins/local/ao_extra_basis/prov_fit_1s.irp.f b/plugins/local/ao_extra_basis/prov_fit_1s.irp.f
index 4d6f07ae..99a18a6b 100644
--- a/plugins/local/ao_extra_basis/prov_fit_1s.irp.f
+++ b/plugins/local/ao_extra_basis/prov_fit_1s.irp.f
@@ -84,20 +84,31 @@ BEGIN_PROVIDER [ integer, list_2p_functions, (n_2p_func_orig)]
  enddo
 END_PROVIDER 
 
+BEGIN_PROVIDER [ integer, extra_fictious_nucl]
+ implicit none
+ extra_fictious_nucl = n_2p_func_tot
+END_PROVIDER 
+
 BEGIN_PROVIDER [ integer, new_nucl_num]
  implicit none
  new_nucl_num = nucl_num + n_2p_func_tot
  print*,'new_nucl_num = ',new_nucl_num
 END_PROVIDER 
 
-BEGIN_PROVIDER [ character*(32), new_nucl_label_1s , (new_nucl_num) ]
+ BEGIN_PROVIDER [ character*(32), new_nucl_label_1s , (new_nucl_num) ]
+&BEGIN_PROVIDER [ integer, list_real_nucl, (nucl_num) ]
+&BEGIN_PROVIDER [ integer, list_fict_nucl, (extra_fictious_nucl) ]
  implicit none
- integer :: i
+ integer :: i,j
  do i = 1, nucl_num 
   new_nucl_label_1s(i) = nucl_label(i)
+  list_real_nucl(i) = i
  enddo
+ j=0
  do i = nucl_num+1,new_nucl_num
+  j+=1
   new_nucl_label_1s(i) = "X"
+  list_fict_nucl(j) = i
  enddo
 END_PROVIDER 
  
@@ -111,6 +122,7 @@ END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, new_nucl_coord_1s_transp, (3,new_nucl_num)]
 &BEGIN_PROVIDER [ double precision, new_nucl_charge_1s, (new_nucl_num)]
+&BEGIN_PROVIDER [ integer, extra_nucl_real_fictious_list_prov, (extra_fictious_nucl)]
  implicit none
  BEGIN_DOC
 ! the real atoms are located in the first nucl_num entries 
@@ -138,11 +150,13 @@ END_PROVIDER
      new_nucl_coord_1s_transp(1:3,k)= nucl_coord_transp(1:3,i)
      new_nucl_coord_1s_transp(good_i,k)+= ao_extra_center
      new_nucl_charge_1s(k) = 0.d0
+     extra_nucl_real_fictious_list_prov(k-nucl_num)=i
      k+=1
      ! one is centered in R_x - d 
      new_nucl_coord_1s_transp(1:3,k)= nucl_coord_transp(1:3,i)
      new_nucl_coord_1s_transp(good_i,k)-= ao_extra_center
      new_nucl_charge_1s(k) = 0.d0
+     extra_nucl_real_fictious_list_prov(k-nucl_num)=i
     enddo
    else if(ao_l(i_ao).gt.1)then
     print*,'WARNING ! Lmax value not implemented yet !'
diff --git a/plugins/local/ao_extra_basis/qp_add_extra_fit_system b/plugins/local/ao_extra_basis/qp_add_extra_fit_system
new file mode 100755
index 00000000..dae95062
--- /dev/null
+++ b/plugins/local/ao_extra_basis/qp_add_extra_fit_system
@@ -0,0 +1,19 @@
+#!/bin/bash 
+# specify the QP folder 
+QP=$QP_ROOT
+dir=${QP}
+# sourcing the quantum_package.rc file
+. ${QP}/quantum_package.rc
+# The main EZFIO folder 
+main=${1%/} 
+# The XYZ file containing the geometry of the additional system you want to add
+extra=${2%.xyz}
+basis_extra=sto-3g
+ezfio_extra=${extra}_${basis_extra}_1s
+echo $ezfio_extra
+qp_fit_1s_basis $extra $basis_extra
+qp set_file $ezfio_extra
+qp run scf | tee ${ezfio_extra}.scf.out
+qp run save_one_e_dm | tee ${ezfio_extra}.one_rdm.out
+qp_copy_extra_basis ${ezfio_extra} $main
+
diff --git a/plugins/local/ao_extra_basis/qp_copy_extra_basis b/plugins/local/ao_extra_basis/qp_copy_extra_basis
index 551bbe1f..cb435e18 100755
--- a/plugins/local/ao_extra_basis/qp_copy_extra_basis
+++ b/plugins/local/ao_extra_basis/qp_copy_extra_basis
@@ -4,7 +4,11 @@ QP=$QP_ROOT
 dir=${QP}
 # sourcing the quantum_package.rc file
 . ${QP}/quantum_package.rc
+# script that copy all data from |AO| basis and nuclei of EZFIO_extra to the ao_extra_basis and extra_nuclei of EZFIO_target
+# use:
+# qp_copy_extra_basis EZFIO_extra EZFIO_target
 EZFIO_extra=${1%/}
+EZFIO_extra=${EZFIO_extra%.xyz}
 EZFIO_target=${2%/}
 
 
@@ -56,4 +60,8 @@ i=primitives_normalized
 newfile=primitives_normalized_extra
 cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
 
+echo "COPYING ALL DATA FROM "$EZFIO_extra"/aux_quantities/ to "${EZFIO_target}"/ao_extra_basis/"
+i=data_one_e_dm_tot_ao.gz
+newfile=ao_extra_one_e_dm.gz
+cp ${EZFIO_extra}/aux_quantities/$i ${EZFIO_target}/ao_extra_basis/$newfile
 
diff --git a/plugins/local/ao_extra_basis/qp_copy_extra_basis_to_usual_basis b/plugins/local/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
index 59be0ee2..ba7e6a71 100755
--- a/plugins/local/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
+++ b/plugins/local/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
@@ -1,5 +1,6 @@
 #!/bin/bash 
 # specify the QP folder 
+# script that copy all data in ao_extra_basis and extra_nuclei and copy it to ao_basis and ao_nuclei
 QP=$QP_ROOT
 dir=${QP}
 # sourcing the quantum_package.rc file
@@ -52,8 +53,5 @@ do
  echo $newfile
  cp ${EZFIO_extra}/ao_extra_basis/$i ${EZFIO_target}/ao_basis/$newfile
 done
-#i=primitives_normalized
-#newfile=primitives_normalized_extra
-#cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
 
 
diff --git a/plugins/local/ao_extra_basis/qp_fit_1s_basis b/plugins/local/ao_extra_basis/qp_fit_1s_basis
new file mode 100755
index 00000000..366e478a
--- /dev/null
+++ b/plugins/local/ao_extra_basis/qp_fit_1s_basis
@@ -0,0 +1,12 @@
+#!/bin/bash 
+## Takes as an argument an xyz file and a basis, and fit the AO basis onto an "s" basis only
+## use: 
+# qp_fit_1s_basis my_xyz_file.xyz basis
+source ~/qp2/quantum_package.rc
+input=${1%.xyz}
+basis=$2
+ezfio_fit=${input}_${basis}_1s
+qp create_ezfio -b $basis $input.xyz -o $ezfio_fit
+# Fitting the original basis on 1s only basis functions
+qp run fit_1s_basis |tee ${ezfio_fit}.fit_1s.out
+qp_copy_extra_basis_to_usual_basis $ezfio_fit $ezfio_fit
diff --git a/plugins/local/ao_extra_basis/transform_basis.irp.f b/plugins/local/ao_extra_basis/transform_basis.irp.f
new file mode 100644
index 00000000..374d2dfc
--- /dev/null
+++ b/plugins/local/ao_extra_basis/transform_basis.irp.f
@@ -0,0 +1,49 @@
+subroutine rotate_nuclei(phi,theta,psi,nucl_centers,n_nucl,nucl_centers_after)
+ implicit none
+ BEGIN_DOC
+ ! routine that rotates a set of nuclei according to three axis corresponding to angles phi, theta, psi.
+ END_DOC
+ double precision, intent(in) :: phi,theta,psi
+ double precision, intent(in) :: nucl_centers(3,n_nucl)
+ integer, intent(in) :: n_nucl
+ double precision, intent(out):: nucl_centers_after(3,n_nucl)
+ double precision :: r_mat(3,3)
+ call r_phi_theta_psi_matrix(phi,theta,psi,r_mat)
+ call get_AB_prod(r_mat,3,3,nucl_centers,n_nucl,nucl_centers_after)
+
+end
+
+
+subroutine r_phi_theta_psi_matrix(phi,theta,psi,r_mat)
+ implicit none
+ BEGIN_DOC
+ ! routine that creates the rotation matrix corresponding to phi,theta,psi 
+ !
+ ! according to conventions in MDFT code
+ END_DOC
+ double precision, intent(in) :: phi,theta,psi
+ double precision, intent(out):: r_mat(3,3)
+ double precision :: ctheta, stheta 
+ double precision :: cphi  , sphi 
+ double precision :: cpsi  , spsi 
+ ctheta = dcos(theta)
+ cphi = dcos(phi)
+ cpsi = dcos(psi)
+
+ stheta = dsin(theta)
+ sphi = dsin(phi)
+ spsi = dsin(psi)
+
+ r_mat(1,1) = ctheta*cphi*cpsi-sphi*spsi
+ r_mat(1,2) = -ctheta*cphi*spsi-sphi*cpsi
+ r_mat(1,3) = stheta*cphi
+ 
+ r_mat(2,1) = ctheta*sphi*cpsi+cphi*spsi
+ r_mat(2,2) = -ctheta*sphi*spsi+cphi*cpsi
+ r_mat(2,3) = stheta*sphi
+
+ r_mat(3,1) = -stheta*cpsi
+ r_mat(3,2) = stheta*spsi
+ r_mat(3,3) = ctheta
+
+end
diff --git a/plugins/local/ao_extra_basis/uninstall b/plugins/local/ao_extra_basis/uninstall
index ff7f8c6a..0e4cd253 100755
--- a/plugins/local/ao_extra_basis/uninstall
+++ b/plugins/local/ao_extra_basis/uninstall
@@ -9,7 +9,7 @@ then
 fi
 
 # list of the scripts to be used by the module 
-scripts_list="qp_copy_extra_basis"
+scripts_list="qp_copy_extra_basis qp_add_extra_fit_system qp_copy_extra_basis_to_usual_basis qp_fit_1s_basis"
 
 # Destroy ONLY the symbolic link for the scripts to be used in the
 # ${QP_ROOT}/scripts/ directory. 
diff --git a/plugins/local/extra_nuclei/EZFIO.cfg b/plugins/local/extra_nuclei/EZFIO.cfg
index e09a7dfd..8c04b4c7 100644
--- a/plugins/local/extra_nuclei/EZFIO.cfg
+++ b/plugins/local/extra_nuclei/EZFIO.cfg
@@ -3,6 +3,12 @@ doc: Number of nuclei
 type: integer
 interface: ezfio, provider
 
+[extra_nucl_pouet]
+doc: Number of nuclei
+type: integer
+interface: ezfio, provider, ocaml
+default:1
+
 [extra_nucl_label]
 doc:  Nuclear labels
 type: character*(32)
@@ -21,3 +27,30 @@ type: double precision
 size: (extra_nuclei.extra_nucl_num,3)
 interface: ezfio
 
+[extra_nucl_real_num]
+doc: Number of real nuclei
+type: integer
+interface: ezfio, provider 
+
+[extra_nucl_fictious_num]
+doc: Number of fictious nuclei
+type: integer
+interface: ezfio, provider 
+
+[extra_nucl_real_fictious_list]
+doc: List of real nuclei to which fictious nuclei are attached to 
+type: integer
+size: (extra_nuclei.extra_nucl_fictious_num)
+interface: ezfio, provider 
+
+[extra_nucl_fictious_list]
+doc: List of fictious nuclei 
+type: integer
+size: (extra_nuclei.extra_nucl_fictious_num)
+interface: ezfio, provider 
+
+[extra_nucl_real_list]
+doc: List of real nuclei
+type: integer
+size: (extra_nuclei.extra_nucl_real_num)
+interface: ezfio, provider 
diff --git a/src/aux_quantities/EZFIO.cfg b/src/aux_quantities/EZFIO.cfg
index 6b4bd0f2..3f6d7783 100644
--- a/src/aux_quantities/EZFIO.cfg
+++ b/src/aux_quantities/EZFIO.cfg
@@ -38,3 +38,10 @@ type: double precision
 size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states)
 
 
+[data_one_e_dm_tot_ao]
+interface: ezfio, provider
+doc: TOTAL (Alpha+Beta) one body density matrix on the |AO| basis computed with the wave function
+type: double precision 
+size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states)
+
+
diff --git a/src/kohn_sham/print_mos.irp.f b/src/kohn_sham/print_mos.irp.f
index 19bb98bc..44718964 100644
--- a/src/kohn_sham/print_mos.irp.f
+++ b/src/kohn_sham/print_mos.irp.f
@@ -1,24 +1,41 @@
 program print_mos
  implicit none
  integer :: i,nx
- double precision :: r(3), xmax, dx, accu
- double precision, allocatable :: mos_array(:)
- double precision:: alpha,envelop,dm_a,dm_b
- allocate(mos_array(mo_num))
- xmax = 5.d0
- nx = 1000
- dx=xmax/dble(nx)
- r = 0.d0
- alpha = 0.5d0
+ print*,ao_kinetic_integrals(1,1)
+ double precision :: xmin,xmax,dx,x,accu,gtt,g,alpha,pi,accu_norm
+ pi = dacos(-1.d0)
+ alpha = 2.D0
+ xmin=0.d0
+ xmax=5.D0
+ nx=10000
+ dx=(xmax-xmin)/dble(nx)
+ x = 0.d0
+ accu = 0.d0
+ accu_norm = 0.d0
  do i = 1, nx
-  call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
-  call give_all_mos_at_r(r,mos_array)
-  accu = mos_array(3)**2+mos_array(4)**2+mos_array(5)**2
-  accu = dsqrt(accu)
-  envelop = (1.d0 - dexp(-alpha * r(3)**2))
-  write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, dm_a+dm_b, envelop
-  r(3) += dx
+  accu += g(x,alpha)*gtt(x,alpha)*x**2 * dx
+  accu_norm += g(x,alpha)**2*x**2 * dx
+  write(33,*)x,g(x,alpha),gtt(x,alpha)
+  x+=dx
  enddo
+ accu=accu * 4.d0 * pi
+ accu_norm *= 4.d0 * pi
+ print*,'accu_norm = ',accu_norm
+ accu*=-0.5D0
+ print*,'accu = ',accu/accu_norm
+end
+
+
+double precision function gtt(x,alpha)
+ implicit none
+ double precision, intent(in) :: x,alpha
+ gtt = dexp(-alpha*x*x) * (4.D0*alpha**2*x**2 - 4.d0 * alpha)
 
 end
 
+double precision function g(x,alpha)
+ implicit none
+ double precision, intent(in) :: x,alpha
+ g = dexp(-alpha*x*x) 
+
+end
diff --git a/src/tools/save_one_e_dm.irp.f b/src/tools/save_one_e_dm.irp.f
index c888f55c..c206343b 100644
--- a/src/tools/save_one_e_dm.irp.f
+++ b/src/tools/save_one_e_dm.irp.f
@@ -27,4 +27,5 @@ subroutine routine_save_one_e_dm
  call ezfio_set_aux_quantities_data_one_e_dm_beta_mo(one_e_dm_mo_beta)
  call ezfio_set_aux_quantities_data_one_e_dm_alpha_ao(one_e_dm_ao_alpha)
  call ezfio_set_aux_quantities_data_one_e_dm_beta_ao(one_e_dm_ao_beta)
+ call ezfio_set_aux_quantities_data_one_e_dm_tot_ao(one_e_dm_ao)
 end