diff --git a/plugins/local/tc_scf/README.md b/plugins/local/tc_scf/README.md
new file mode 100644
index 00000000..45f5c60f
--- /dev/null
+++ b/plugins/local/tc_scf/README.md
@@ -0,0 +1,2 @@
+Transcorrelated SCF
+===================
diff --git a/src/ao_extra_basis/EZFIO.cfg b/src/ao_extra_basis/EZFIO.cfg
new file mode 100644
index 00000000..8841c194
--- /dev/null
+++ b/src/ao_extra_basis/EZFIO.cfg
@@ -0,0 +1,104 @@
+[ao_extra_basis]
+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|
+interface: ezfio, provider
+
+[ao_extra_prim_num]
+type: integer
+doc: Number of primitives per |ao_extra|
+size: (ao_extra_basis.ao_extra_num)
+interface: ezfio, provider
+
+[ao_extra_prim_num_max]
+type: integer
+doc: Maximum number of primitives
+default: =maxval(ao_extra_basis.ao_extra_prim_num)
+interface: ezfio
+
+[ao_extra_nucl]
+type: integer
+doc: Index of the nucleus on which the |ao_extra| is centered
+size: (ao_extra_basis.ao_extra_num)
+interface: ezfio, provider
+
+[ao_extra_power]
+type: integer
+doc: Powers of x, y and z for each |ao_extra|
+size: (ao_extra_basis.ao_extra_num,3)
+interface: ezfio, provider
+
+[ao_extra_coef]
+type: double precision
+doc: Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** ao_extras.
+size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
+interface: ezfio, provider
+
+[ao_extra_expo]
+type: double precision
+doc: Exponents for each primitive of each |ao_extra|
+size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
+interface: ezfio, provider
+
+[ao_extra_md5]
+type: character*(32)
+doc: MD5 key, specific of the |ao_extra| basis
+interface: ezfio, provider
+
+[ao_extra_cartesian]
+type: logical
+doc: If |true|, use |ao_extras| in Cartesian coordinates (6d,10f,...)
+interface: ezfio, provider
+default: false
+
+[ao_extra_normalized]
+type: logical
+doc: Use normalized basis functions
+interface: ezfio, provider
+default: true
+
+[primitives_normalized_extra]
+type: logical
+doc: Use normalized primitive functions
+interface: ezfio, provider
+default: true
+
+[ao_extra_expo_im]
+type: double precision
+doc: imag part for Exponents for each primitive of each cGTOs |ao_extra|
+size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
+interface: ezfio, provider
+
+[ao_extra_expo_pw]
+type: double precision
+doc: plane wave part 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_expo_phase]
+type: double precision
+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/src/ao_extra_basis/LiH.xyz b/src/ao_extra_basis/LiH.xyz
new file mode 100644
index 00000000..03431803
--- /dev/null
+++ b/src/ao_extra_basis/LiH.xyz
@@ -0,0 +1,4 @@
+2
+
+H 0. 0. 0.
+Li 0. 0. 1.0
diff --git a/src/ao_extra_basis/NEED b/src/ao_extra_basis/NEED
new file mode 100644
index 00000000..c4e8c3cf
--- /dev/null
+++ b/src/ao_extra_basis/NEED
@@ -0,0 +1,3 @@
+extra_nuclei
+basis
+ao_basis
diff --git a/src/ao_extra_basis/README.rst b/src/ao_extra_basis/README.rst
new file mode 100644
index 00000000..f60d71c0
--- /dev/null
+++ b/src/ao_extra_basis/README.rst
@@ -0,0 +1,15 @@
+===========
+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/src/ao_extra_basis/aos.irp.f b/src/ao_extra_basis/aos.irp.f
new file mode 100644
index 00000000..56d6fb04
--- /dev/null
+++ b/src/ao_extra_basis/aos.irp.f
@@ -0,0 +1,325 @@
+BEGIN_PROVIDER [ integer, ao_extra_prim_num_max ]
+ implicit none
+ BEGIN_DOC
+ ! Max number of primitives.
+ END_DOC
+ ao_extra_prim_num_max = maxval(ao_extra_prim_num)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_extra_shell, (ao_extra_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the ao_extra corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)+1
+   do j=1,(n*(n+1))/2
+     k = k+1
+     ao_extra_shell(k) = i
+   enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_extra_first_of_shell, (shell_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the ao_extra corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=1
+ do i=1,shell_num
+   ao_extra_first_of_shell(i) = k
+   n = shell_ang_mom(i)+1
+   k = k+(n*(n+1))/2
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized, (ao_extra_num,ao_extra_prim_num_max) ]
+&BEGIN_PROVIDER [ double precision, ao_extra_coef_normalization_factor, (ao_extra_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Coefficients including the |ao_extra| normalization
+  END_DOC
+  double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+  integer                        :: l, powA(3), nz
+  integer                        :: i,j,k
+  nz=100
+  C_A(1) = 0.d0
+  C_A(2) = 0.d0
+  C_A(3) = 0.d0
+  ao_extra_coef_normalized = 0.d0
+
+  do i=1,ao_extra_num
+
+!    powA(1) = ao_extra_power(i,1) +  ao_extra_power(i,2) +  ao_extra_power(i,3)
+!    powA(2) = 0
+!    powA(3) = 0
+    powA(1) = ao_extra_power(i,1)
+    powA(2) = ao_extra_power(i,2)
+    powA(3) = ao_extra_power(i,3)
+
+    ! Normalization of the primitives
+    if (primitives_normalized_extra) then
+      do j=1,ao_extra_prim_num(i)
+        call overlap_gaussian_xyz(C_A,C_A,ao_extra_expo(i,j),ao_extra_expo(i,j), &
+           powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+        ao_extra_coef_normalized(i,j) = ao_extra_coef(i,j)/dsqrt(norm)
+      enddo
+    else
+      do j=1,ao_extra_prim_num(i)
+        ao_extra_coef_normalized(i,j) = ao_extra_coef(i,j)
+      enddo
+    endif
+
+    ! Normalization of the contracted basis functions
+    if (ao_extra_normalized) then
+      norm = 0.d0
+      do j=1,ao_extra_prim_num(i)
+        do k=1,ao_extra_prim_num(i)
+          call overlap_gaussian_xyz(C_A,C_A,ao_extra_expo(i,j),ao_extra_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+          norm = norm+c*ao_extra_coef_normalized(i,j)*ao_extra_coef_normalized(i,k)
+        enddo
+      enddo
+      ao_extra_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
+    else
+      ao_extra_coef_normalization_factor(i) = 1.d0
+    endif
+  enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered, (ao_extra_num,ao_extra_prim_num_max) ]
+&BEGIN_PROVIDER [ double precision, ao_extra_expo_ordered, (ao_extra_num,ao_extra_prim_num_max) ]
+  implicit none
+  BEGIN_DOC
+  ! Sorted primitives to accelerate 4 index |MO| transformation
+  END_DOC
+
+  integer                        :: iorder(ao_extra_prim_num_max)
+  double precision               :: d(ao_extra_prim_num_max,2)
+  integer                        :: i,j
+  do i=1,ao_extra_num
+    do j=1,ao_extra_prim_num(i)
+      iorder(j) = j
+      d(j,1) = ao_extra_expo(i,j)
+      d(j,2) = ao_extra_coef_normalized(i,j)
+    enddo
+    call dsort(d(1,1),iorder,ao_extra_prim_num(i))
+    call dset_order(d(1,2),iorder,ao_extra_prim_num(i))
+    do j=1,ao_extra_prim_num(i)
+      ao_extra_expo_ordered(i,j) = d(j,1)
+      ao_extra_coef_normalized_ordered(i,j) = d(j,2)
+    enddo
+  enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered_transp, (ao_extra_prim_num_max,ao_extra_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Transposed :c:data:`ao_extra_coef_normalized_ordered`
+  END_DOC
+  integer                        :: i,j
+  do j=1, ao_extra_num
+    do i=1, ao_extra_prim_num_max
+      ao_extra_coef_normalized_ordered_transp(i,j) = ao_extra_coef_normalized_ordered(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_extra_expo_ordered_transp, (ao_extra_prim_num_max,ao_extra_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Transposed :c:data:`ao_extra_expo_ordered`
+  END_DOC
+  integer                        :: i,j
+  do j=1, ao_extra_num
+    do i=1, ao_extra_prim_num_max
+      ao_extra_expo_ordered_transp(i,j) = ao_extra_expo_ordered(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, ao_extra_l, (ao_extra_num) ]
+&BEGIN_PROVIDER [ integer, ao_extra_l_max  ]
+&BEGIN_PROVIDER [ character*(128), ao_extra_l_char, (ao_extra_num) ]
+ implicit none
+ BEGIN_DOC
+! :math:`l` value of the |ao_extra|: :math`a+b+c` in :math:`x^a y^b z^c`
+ END_DOC
+ integer :: i
+ do i=1,ao_extra_num
+   ao_extra_l(i) = ao_extra_power(i,1) + ao_extra_power(i,2) + ao_extra_power(i,3)
+   ao_extra_l_char(i) = l_to_character(ao_extra_l(i))
+ enddo
+ ao_extra_l_max = maxval(ao_extra_l)
+END_PROVIDER
+
+integer function ao_extra_power_index(nx,ny,nz)
+  implicit none
+  integer, intent(in)            :: nx, ny, nz
+  BEGIN_DOC
+  ! Unique index given to a triplet of powers:
+  !
+  ! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1`
+  END_DOC
+  integer                        :: l
+  l = nx + ny + nz
+  ao_extra_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
+end
+
+
+ BEGIN_PROVIDER [ integer, Nucl_N_ao_extras, (extra_nucl_num)]
+&BEGIN_PROVIDER [ integer, N_ao_extras_max ]
+ implicit none
+ BEGIN_DOC
+ ! Number of |ao_extras| per atom
+ END_DOC
+ integer :: i
+ Nucl_N_ao_extras = 0
+ do i = 1, ao_extra_num
+  Nucl_N_ao_extras(ao_extra_nucl(i)) +=1
+ enddo
+ N_ao_extras_max = maxval(Nucl_N_ao_extras)
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, Nucl_ao_extras, (extra_nucl_num,N_ao_extras_max)]
+ implicit none
+ BEGIN_DOC
+ ! List of |ao_extras| centered on each atom
+ END_DOC
+ integer :: i
+ integer, allocatable :: nucl_tmp(:)
+ allocate(nucl_tmp(nucl_num))
+ nucl_tmp = 0
+ Nucl_ao_extras = 0
+ do i = 1, ao_extra_num
+  nucl_tmp(ao_extra_nucl(i))+=1
+  Nucl_ao_extras(ao_extra_nucl(i),nucl_tmp(ao_extra_nucl(i))) = i
+ enddo
+ deallocate(nucl_tmp)
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, Nucl_list_shell_ao_extras, (extra_nucl_num,N_ao_extras_max)]
+&BEGIN_PROVIDER [ integer, Nucl_num_shell_ao_extras, (nucl_num)]
+ implicit none
+ integer :: i,j,k
+ BEGIN_DOC
+ ! Index of the shell type |ao_extras| and of the corresponding |ao_extras|
+ ! By convention, for p,d,f and g |ao_extras|, we take the index
+ ! of the |ao_extra| with the the corresponding power in the x axis
+ END_DOC
+ do i = 1, extra_nucl_num
+  Nucl_num_shell_ao_extras(i) = 0
+  do j = 1, Nucl_N_ao_extras(i)
+    if (ao_extra_power(Nucl_ao_extras(i,j),1) == ao_extra_l(Nucl_ao_extras(i,j))) then
+     Nucl_num_shell_ao_extras(i)+=1
+     Nucl_list_shell_ao_extras(i,Nucl_num_shell_ao_extras(i))=Nucl_ao_extras(i,j)
+    endif
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ character*(4), ao_extra_l_char_space, (ao_extra_num) ]
+ implicit none
+ BEGIN_DOC
+! Converts an l value to a string
+ END_DOC
+ integer :: i
+ character*(4) :: give_ao_extra_character_space
+ do i=1,ao_extra_num
+
+  if(ao_extra_l(i)==0)then
+  ! S type ao_extra
+   give_ao_extra_character_space = 'S   '
+  elseif(ao_extra_l(i) == 1)then
+  ! P type ao_extra
+   if(ao_extra_power(i,1)==1)then
+    give_ao_extra_character_space = 'X   '
+   elseif(ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'Y   '
+   else
+    give_ao_extra_character_space = 'Z   '
+   endif
+  elseif(ao_extra_l(i) == 2)then
+  ! D type ao_extra
+   if(ao_extra_power(i,1)==2)then
+    give_ao_extra_character_space = 'XX  '
+   elseif(ao_extra_power(i,2) == 2)then
+    give_ao_extra_character_space = 'YY  '
+   elseif(ao_extra_power(i,3) == 2)then
+    give_ao_extra_character_space = 'ZZ  '
+   elseif(ao_extra_power(i,1) == 1 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'XY  '
+   elseif(ao_extra_power(i,1) == 1 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'XZ  '
+   else
+    give_ao_extra_character_space = 'YZ  '
+   endif
+  elseif(ao_extra_l(i) == 3)then
+  ! F type ao_extra
+   if(ao_extra_power(i,1)==3)then
+    give_ao_extra_character_space = 'XXX '
+   elseif(ao_extra_power(i,2) == 3)then
+    give_ao_extra_character_space = 'YYY '
+   elseif(ao_extra_power(i,3) == 3)then
+    give_ao_extra_character_space = 'ZZZ '
+   elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'XXY '
+   elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'XXZ '
+   elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,1) == 1)then
+    give_ao_extra_character_space = 'YYX '
+   elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'YYZ '
+   elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,1) == 1)then
+    give_ao_extra_character_space = 'ZZX '
+   elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'ZZY '
+   elseif(ao_extra_power(i,3) == 1 .and. ao_extra_power(i,2) == 1 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'XYZ '
+   endif
+  elseif(ao_extra_l(i) == 4)then
+  ! G type ao_extra
+   if(ao_extra_power(i,1)==4)then
+    give_ao_extra_character_space = 'XXXX'
+   elseif(ao_extra_power(i,2) == 4)then
+    give_ao_extra_character_space = 'YYYY'
+   elseif(ao_extra_power(i,3) == 4)then
+    give_ao_extra_character_space = 'ZZZZ'
+   elseif(ao_extra_power(i,1) == 3 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'XXXY'
+   elseif(ao_extra_power(i,1) == 3 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'XXXZ'
+   elseif(ao_extra_power(i,2) == 3 .and. ao_extra_power(i,1) == 1)then
+    give_ao_extra_character_space = 'YYYX'
+   elseif(ao_extra_power(i,2) == 3 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'YYYZ'
+   elseif(ao_extra_power(i,3) == 3 .and. ao_extra_power(i,1) == 1)then
+    give_ao_extra_character_space = 'ZZZX'
+   elseif(ao_extra_power(i,3) == 3 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'ZZZY'
+   elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 2)then
+    give_ao_extra_character_space = 'XXYY'
+   elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,3) == 2)then
+    give_ao_extra_character_space = 'YYZZ'
+   elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 1 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'XXYZ'
+   elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,1) == 1 .and. ao_extra_power(i,3) == 1)then
+    give_ao_extra_character_space = 'YYXZ'
+   elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,1) == 1 .and. ao_extra_power(i,2) == 1)then
+    give_ao_extra_character_space = 'ZZXY'
+   endif
+  endif
+   ao_extra_l_char_space(i) = give_ao_extra_character_space
+ enddo
+END_PROVIDER
diff --git a/src/ao_extra_basis/aos_transp.irp.f b/src/ao_extra_basis/aos_transp.irp.f
new file mode 100644
index 00000000..ed34835b
--- /dev/null
+++ b/src/ao_extra_basis/aos_transp.irp.f
@@ -0,0 +1,68 @@
+
+! ---
+
+BEGIN_PROVIDER [ integer, Nucl_ao_extras_transposed, (N_ao_extras_max,nucl_num)]
+
+  BEGIN_DOC
+  ! List of ao_extras attached on each atom
+  END_DOC
+ 
+  implicit none
+  integer              :: i
+  integer, allocatable :: nucl_tmp(:)
+
+  allocate(nucl_tmp(nucl_num))
+  nucl_tmp = 0
+  do i = 1, ao_extra_num
+    nucl_tmp(ao_extra_nucl(i)) += 1
+    Nucl_ao_extras_transposed(nucl_tmp(ao_extra_nucl(i)),ao_extra_nucl(i)) = i
+  enddo
+  deallocate(nucl_tmp)
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, ao_extra_expo_ordered_transp_per_nucl, (ao_extra_prim_num_max,N_ao_extras_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_ao_extras(i)
+   k = Nucl_ao_extras_transposed(j,i)
+   do l = 1, ao_extra_prim_num(k)
+    ao_extra_expo_ordered_transp_per_nucl(l,j,i) = ao_extra_expo_ordered_transp(l,k)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, ao_extra_power_ordered_transp_per_nucl, (3,N_ao_extras_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_ao_extras(i)
+   k = Nucl_ao_extras_transposed(j,i)
+   do l = 1, 3
+    ao_extra_power_ordered_transp_per_nucl(l,j,i) = ao_extra_power(k,l)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered_transp_per_nucl, (ao_extra_prim_num_max,N_ao_extras_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_ao_extras(i)
+   k = Nucl_ao_extras_transposed(j,i)
+   do l = 1, ao_extra_prim_num(k)
+    ao_extra_coef_normalized_ordered_transp_per_nucl(l,j,i) = ao_extra_coef_normalized_ordered_transp(l,k)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
diff --git a/src/ao_extra_basis/density_extra.irp.f b/src/ao_extra_basis/density_extra.irp.f
new file mode 100644
index 00000000..5d6492e3
--- /dev/null
+++ b/src/ao_extra_basis/density_extra.irp.f
@@ -0,0 +1,87 @@
+BEGIN_PROVIDER [ double precision, effective_ao_extra_dm, (ao_extra_num, ao_extra_num)]
+ implicit none
+ BEGIN_DOC
+ ! effective density matrix : rho_pq x N_p x N_q x E_pq x (pi/gamma_pq)^3/2
+ !
+ ! where rho_pq is the usual density matrix 
+ !
+ ! N_p and N_q are the normalization factors associated with the two Gaussians 
+ !
+ ! E_pq is the prefactor resulting from the Gaussian product 
+ !
+ ! gamma_pq = gamma_p + gamm_q
+ END_DOC
+ integer :: i,j
+ do i = 1, ao_extra_num
+  do j = 1, ao_extra_num
+   effective_ao_extra_dm(j,i) = ao_extra_one_e_dm(j,i) * ao_extra_coef_normalized(j,1) * ao_extra_coef_normalized(i,1) & 
+       * inv_pi_gamma_pq_3_2_ao_extra(j,i) * E_pq_ao_extra(j,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
+&BEGIN_PROVIDER [ double precision, inv_gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
+&BEGIN_PROVIDER [ double precision, inv_pi_gamma_pq_3_2_ao_extra, (ao_extra_num,ao_extra_num)]
+&BEGIN_PROVIDER [ double precision, sqrt_gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
+ implicit none
+ BEGIN_DOC
+ ! gamma_pq_ao_extra = gamma_p + gamma_q  
+ !
+ ! inv_gamma_pq_ao_extra = 1/(gamma_p + gamma_q)
+ !
+ ! sqrt_gamma_pq_ao_extra = sqrt(gamma_p + gamma_q)
+ !
+ ! WARNING :: VALID ONLY IN THE CASE OF A PURELY 1S BASIS 
+ END_DOC
+ include 'constants.include.F'
+ integer :: i,j
+ do i = 1, ao_extra_num
+  do j = 1, ao_extra_num
+   gamma_pq_ao_extra(j,i) = ao_extra_expo(j,1) + ao_extra_expo(i,1)
+   inv_gamma_pq_ao_extra(j,i) = 1.d0/gamma_pq_ao_extra(j,i)
+   sqrt_gamma_pq_ao_extra(j,i) = dsqrt(gamma_pq_ao_extra(j,i))
+   inv_pi_gamma_pq_3_2_ao_extra(j,i) = (pi * inv_gamma_pq_ao_extra(j,i))**(1.5d0)
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, E_pq_ao_extra, (ao_extra_num,ao_extra_num)]
+ implicit none
+ BEGIN_DOC
+ ! E_pq_ao_extra = exp(-alpha_p alpha_q/gamma_pq (Q_p - Q_q)^2)
+ END_DOC 
+ integer :: i,j
+ do i = 1, ao_extra_num
+  do j = 1, ao_extra_num
+   E_pq_ao_extra(j,i) = dexp(-ao_extra_center_1s_dist(j,i)**2 * ao_extra_expo(j,1)*ao_extra_expo(i,1)*inv_gamma_pq_ao_extra(j,i))
+  enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, ao_extra_center_1s, (3,ao_extra_num)]
+ implicit none
+ BEGIN_DOC
+ ! Original position of each extra AO 
+ END_DOC
+ integer :: i,i_nucl
+ do i = 1, ao_extra_num
+  i_nucl= ao_extra_nucl(i)
+  ao_extra_center_1s(1:3,i) = extra_nucl_coord(i_nucl,1:3)
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, ao_extra_center_1s_dist, (ao_extra_num,ao_extra_num)]
+ implicit none
+ integer :: i,j
+ do i = 1, ao_extra_num
+  do j = 1, ao_extra_num
+   ao_extra_center_1s_dist(j,i) = (ao_extra_center_1s(1,j) - ao_extra_center_1s(1,i))**2 &
+                                + (ao_extra_center_1s(2,j) - ao_extra_center_1s(2,i))**2 &
+                                + (ao_extra_center_1s(3,j) - ao_extra_center_1s(3,i))**2
+    ao_extra_center_1s_dist(j,i)=dsqrt(ao_extra_center_1s_dist(j,i))
+  enddo
+ enddo
+END_PROVIDER 
diff --git a/src/ao_extra_basis/dimensions_integrals.irp.f b/src/ao_extra_basis/dimensions_integrals.irp.f
new file mode 100644
index 00000000..80dcda18
--- /dev/null
+++ b/src/ao_extra_basis/dimensions_integrals.irp.f
@@ -0,0 +1,19 @@
+ BEGIN_PROVIDER [ integer, n_pt_max_extra_basis_integrals ]
+&BEGIN_PROVIDER [ integer, n_pt_max_extra_basis_i_x]
+ implicit none
+ BEGIN_DOC
+! Number of points used in the numerical integrations.
+ END_DOC
+ integer :: n_pt_sup
+ integer :: prim_power_l_max
+ include 'utils/constants.include.F'
+ prim_power_l_max = maxval(ao_extra_power)
+ n_pt_max_extra_basis_integrals = 24 * prim_power_l_max + 4
+ n_pt_max_extra_basis_i_x = 8 * prim_power_l_max
+ ASSERT (n_pt_max_extra_basis_i_x-1 <= max_dim)
+ if (n_pt_max_extra_basis_i_x-1 > max_dim) then
+   print *, 'Increase max_dim in utils/constants.include.F to ', n_pt_max_extra_basis_i_x-1
+   stop 1
+ endif
+END_PROVIDER
+
diff --git a/src/ao_extra_basis/extra_basis.irp.f b/src/ao_extra_basis/extra_basis.irp.f
new file mode 100644
index 00000000..9773af2f
--- /dev/null
+++ b/src/ao_extra_basis/extra_basis.irp.f
@@ -0,0 +1,16 @@
+program extra_basis
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+ integer :: i
+ print*,'extra_nucl_num = ',extra_nucl_num
+ do i = 1, extra_nucl_num
+  print*,'i = ',i
+  print*,'extra_nucl_label  = ',extra_nucl_label(i)
+  print*,'extra_nucl_charge = ',extra_nucl_charge(i)
+  print*,'extra_nucl_coord  = '
+  print*,extra_nucl_coord(i,1:3)
+ enddo
+ print*,ao_extra_num
+end
diff --git a/src/ao_extra_basis/fit_1s_basis.irp.f b/src/ao_extra_basis/fit_1s_basis.irp.f
new file mode 100644
index 00000000..ef09d5b2
--- /dev/null
+++ b/src/ao_extra_basis/fit_1s_basis.irp.f
@@ -0,0 +1,41 @@
+program fit_1s_basis
+ implicit none
+ provide lmax_too_big
+ integer :: i,j
+ print*,'////////////////////////////////////////////////////'
+ print*,'////////////////////////////////////////////////////'
+ print*,'Fitting the original basis set on uncontracted s only functions '
+ print*,'WARNING :: works for now with only P functions at most !!'
+ print*,'WARNING :: otherwise it will stop '
+ 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)
+ call ezfio_set_ao_extra_basis_ao_extra_expo(new_ao_expo_1s)
+ call ezfio_set_ao_extra_basis_ao_extra_power(new_ao_power_1s)
+end
+
diff --git a/src/ao_extra_basis/h2o.xyz b/src/ao_extra_basis/h2o.xyz
new file mode 100644
index 00000000..d3928214
--- /dev/null
+++ b/src/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/src/ao_extra_basis/install b/src/ao_extra_basis/install
new file mode 100755
index 00000000..e7e668e1
--- /dev/null
+++ b/src/ao_extra_basis/install
@@ -0,0 +1,23 @@
+#!/bin/bash
+                                                                                                                                                        
+# Check if the QP_ROOT environment variable is set. 
+if [[ -z ${QP_ROOT} ]]
+then
+  print "The QP_ROOT environment variable is not set."
+  print "Please reload the quantum_package.rc file."
+  exit -1
+fi
+
+# Get the absolute path of the current directory.
+currdir=${PWD}
+
+# list of the scripts to be used by the module 
+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. 
+ for i in $scripts_list
+ do
+  ln --symbolic ${currdir}/$i ${QP_ROOT}/scripts/
+ done
+
diff --git a/src/ao_extra_basis/prov_fit_1s.irp.f b/src/ao_extra_basis/prov_fit_1s.irp.f
new file mode 100644
index 00000000..99a18a6b
--- /dev/null
+++ b/src/ao_extra_basis/prov_fit_1s.irp.f
@@ -0,0 +1,266 @@
+BEGIN_PROVIDER [ double precision, ao_extra_center]
+ implicit none
+ ao_extra_center = 0.01d0
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ integer, n_func_tot]
+ implicit none
+ BEGIN_DOC
+ ! n_func_tot :: total number of functions in the fitted basis set
+ ! 
+ ! returned in an uncontracted way
+ END_DOC
+ integer :: i,prefact
+ n_func_tot = 0
+   print*,'n_func_tot '
+ do i = 1, ao_num
+   if(ao_l(i) == 0)then
+    prefact = 1 ! s functions 
+   else
+    ! p functions are fitted with 2 functions
+    ! d functions are fitted with 4 functions etc ...
+    prefact=2*ao_l(i) 
+   endif
+   n_func_tot += prefact * ao_prim_num(i) 
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ integer, n_prim_tot_orig]
+ implicit none
+ integer :: i
+ n_prim_tot_orig = 0
+ do i = 1, ao_num
+  n_prim_tot_orig += ao_prim_num(i) 
+ enddo
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ logical, lmax_too_big]
+ implicit none
+ if (ao_l_max.gt.1)then
+  lmax_too_big = .True.
+ else
+  lmax_too_big = .False.
+ endif
+ if(lmax_too_big)then
+  print*,'STOPPING !! lmax is larger than 1 !'
+  print*,'Cannot yet fit with 1s functions ...'
+  stop
+ endif
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ integer, n_2p_func_orig]
+&BEGIN_PROVIDER [ integer, n_2p_func_tot]
+ implicit none
+ integer :: i
+ BEGIN_DOC
+ ! n_2p_func_orig :: number of 2p functions in the original basis
+ !
+ ! n_2p_func_tot :: total number of p functions in the fitted basis
+ END_DOC
+ n_2p_func_orig= 0
+ n_2p_func_tot = 0
+ do i = 1, ao_num
+  if(ao_l(i)==1)then
+   n_2p_func_orig+= 1
+   n_2p_func_tot += ao_prim_num(i) * 2
+  endif
+ enddo
+ print*,'n_2p_func_tot = ',n_2p_func_tot
+END_PROVIDER 
+
+BEGIN_PROVIDER [ integer, list_2p_functions, (n_2p_func_orig)]
+ implicit none
+ BEGIN_DOC
+ ! list of 2p functions in the original basis
+ END_DOC
+ integer :: i,j
+ j=0
+ do i = 1, ao_num
+  if(ao_l(i)==1)then
+   j+=1
+   list_2p_functions(j) = i
+  endif
+ 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 [ integer, list_real_nucl, (nucl_num) ]
+&BEGIN_PROVIDER [ integer, list_fict_nucl, (extra_fictious_nucl) ]
+ implicit none
+ 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 
+ 
+ BEGIN_PROVIDER [ double precision,  new_nucl_coord_1s, (new_nucl_num,3)]
+ implicit none
+ integer :: i,j
+ do i = 1, new_nucl_num
+  new_nucl_coord_1s(i,1:3) = new_nucl_coord_1s_transp(1:3,i)
+ enddo
+ 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 
+! 
+! then the fictious atoms are located after
+ END_DOC
+ integer :: i,ii,j,i_ao,k,n_ao
+ integer :: return_xyz_int,power(3),good_i
+ new_nucl_charge_1s = 0.d0
+ do i = 1, nucl_num
+  new_nucl_coord_1s_transp(1:3,i) = nucl_coord_transp(1:3,i)
+  new_nucl_charge_1s(i) = nucl_charge(i)
+ enddo
+ k = nucl_num
+ do i = 1, nucl_num
+  do ii = 1, Nucl_N_Aos(i)
+   i_ao = nucl_aos_transposed(ii,i)
+   if(ao_l(i_ao)==1)then
+    ! split the function into 2 s functions 
+    ! one is centered in R_x + d 
+    power(1:3) = ao_power(i_ao,1:3)
+    good_i = return_xyz_int(power)
+    do j = 1, ao_prim_num(i_ao)
+     k+=1
+     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 !'
+    print*,'stopping ...'
+    stop
+   endif
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ integer, new_n_AOs_max]
+ implicit none
+ new_n_AOs_max = ao_prim_num_max * n_AOs_max
+ 
+ END_PROVIDER 
+
+
+ BEGIN_PROVIDER [ integer, new_Nucl_N_Aos, (new_nucl_num)]
+&BEGIN_PROVIDER [ integer, new_nucl_aos_transposed, (new_n_AOs_max,new_nucl_num) ]
+&BEGIN_PROVIDER [ double precision, new_ao_expo_1s , (n_func_tot) ]
+&BEGIN_PROVIDER [ integer, new_ao_nucl_1s, (n_func_tot)]
+ implicit none
+ integer :: i,j,ii,i_ao,n_func,n_func_total,n_nucl
+ double precision :: coef
+ n_func_total = 0
+ do i = 1, nucl_num
+  n_func = 0
+  do ii = 1, Nucl_N_Aos(i)
+   i_ao = nucl_aos_transposed(ii,i)
+   if(ao_l(i_ao)==0)then
+    do j = 1, ao_prim_num(i_ao)
+     coef= ao_expo(i_ao,j)
+     n_func_total += 1
+     n_func +=1
+     new_nucl_aos_transposed(n_func,i) = n_func_total
+     new_ao_expo_1s(n_func_total) = coef
+     new_ao_nucl_1s(n_func_total) = i
+    enddo
+   endif
+  enddo
+  new_Nucl_N_Aos(i) = n_func
+ enddo
+ n_nucl=nucl_num
+ do i = 1, nucl_num
+  do ii = 1, Nucl_N_Aos(i)
+   i_ao = nucl_aos_transposed(ii,i)
+   if(ao_l(i_ao)==1)then
+    do j = 1, ao_prim_num(i_ao)
+     coef= ao_expo(i_ao,j)
+     n_func_total+=1
+     n_nucl +=1
+     new_nucl_aos_transposed(1,n_nucl) = n_func_total
+     new_ao_expo_1s(n_func_total) = coef
+     new_Nucl_N_Aos(n_nucl)=1
+     new_ao_nucl_1s(n_func_total) = n_nucl
+
+     n_func_total+=1
+     n_nucl +=1
+     new_nucl_aos_transposed(1,n_nucl) = n_func_total
+     new_ao_expo_1s(n_func_total) = coef
+     new_Nucl_N_Aos(n_nucl)=1
+     new_ao_nucl_1s(n_func_total) = n_nucl
+    enddo
+   endif
+  enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, new_ao_coef_1s , (n_func_tot) ]
+  implicit none
+ integer :: i
+  BEGIN_DOC
+! Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** AOs.
+  END_DOC
+  do i = 1, n_func_tot
+   new_ao_coef_1s(i) = 1.d0
+  enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ integer, new_ao_prim_num_1s, (n_func_tot)]
+ implicit none
+ integer :: i
+ do i = 1, n_func_tot
+  new_ao_prim_num_1s(i) = 1
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, new_ao_power_1s, (n_func_tot,3)]
+ implicit none
+ new_ao_power_1s = 0
+END_PROVIDER 
+
+integer function return_xyz_int(power)
+ implicit none
+ integer,intent(in) :: power(3)
+ if(power(1) == 1 .and. power(2) ==0 .and. power(3) ==0)then
+  return_xyz_int = 1
+ else if (power(2) == 1 .and. power(1) ==0 .and. power(3) ==0)then
+  return_xyz_int = 2
+ else if (power(3) == 1 .and. power(1) ==0 .and. power(1) ==0)then
+  return_xyz_int = 3
+ else 
+  return_xyz_int = -1000
+ endif
+end
diff --git a/src/ao_extra_basis/qp_add_extra_fit_system b/src/ao_extra_basis/qp_add_extra_fit_system
new file mode 100755
index 00000000..dae95062
--- /dev/null
+++ b/src/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/src/ao_extra_basis/qp_copy_extra_basis b/src/ao_extra_basis/qp_copy_extra_basis
new file mode 100755
index 00000000..cb435e18
--- /dev/null
+++ b/src/ao_extra_basis/qp_copy_extra_basis
@@ -0,0 +1,67 @@
+#!/bin/bash 
+# specify the QP folder 
+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%/}
+
+
+
+echo "********** SCRIPT TO COPY DATA FROM EZFIO TO ANOTHER *********"
+echo "Extracting data from "$EZFIO_extra 
+echo "and Copying data to  "$EZFIO_target
+
+### COPYING ALL DATA FROM $EZFIO_extra/nuclei/ to ${EZFIO_target}/extra_nuclei/
+echo "COPYING ALL DATA FROM "$EZFIO_extra"/nuclei/ to "${EZFIO_target}"/extra_nuclei/"
+direxists=false
+if [ -d ${EZFIO_target}/extra_nuclei/ ] ; then
+ direxists=true
+ echo "The directory extra_nuclei exists"
+ else
+ echo "Creating the directory extra_nuclei "
+ direxists=false
+ mkdir ${EZFIO_target}/extra_nuclei/
+fi
+data=`\ls $EZFIO_extra/nuclei/`
+for i in $data
+do
+ echo $i
+ newfile=`echo $i | sed 's/nucl/extra_nucl/g' `
+ echo $newfile
+ cp ${EZFIO_extra}/nuclei/$i ${EZFIO_target}/extra_nuclei/$newfile
+done
+
+### COPYING ALL DATA FROM $EZFIO_extra/ao_basis/ to ${EZFIO_target}/ao_extra_basis/
+direxists=false
+if [ -d ${EZFIO_target}/ao_extra_basis/ ] ; then
+ direxists=true
+ echo "The directory exists ao_extra_basis"
+ else
+ echo "Creating the directory ao_extra_basis"
+ direxists=false
+ mkdir ${EZFIO_target}/ao_extra_basis/
+fi
+echo "COPYING ALL DATA FROM "$EZFIO_extra"/ao_basis/ to "${EZFIO_target}"/ao_extra_basis/"
+data=`\ls $EZFIO_extra/ao_basis/`
+for i in $data
+do
+ echo $i
+ newfile=`echo $i | sed 's/ao/ao_extra/g' `
+ echo $newfile
+ cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
+done
+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/src/ao_extra_basis/qp_copy_extra_basis_to_usual_basis b/src/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
new file mode 100755
index 00000000..ba7e6a71
--- /dev/null
+++ b/src/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
@@ -0,0 +1,57 @@
+#!/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
+. ${QP}/quantum_package.rc
+EZFIO_extra=${1%/}
+EZFIO_target=${2%/}
+
+
+
+echo "********** SCRIPT TO COPY DATA FROM EZFIO TO ANOTHER *********"
+echo "Extracting data from "$EZFIO_extra 
+echo "and Copying data to  "$EZFIO_target
+
+### COPYING ALL DATA FROM $EZFIO_extra/nuclei/ to ${EZFIO_target}/extra_nuclei/
+echo "COPYING ALL DATA FROM "$EZFIO_extra"/extra_nuclei/ to "${EZFIO_target}"/nuclei/"
+direxists=false
+if [ -d ${EZFIO_target}/extra_nuclei/ ] ; then
+ direxists=true
+ echo "The directory extra_nuclei exists"
+ else
+ echo "PB !!"
+ direxists=false
+ mkdir ${EZFIO_target}/extra_nuclei/
+fi
+data=`\ls $EZFIO_extra/extra_nuclei/`
+for i in $data
+do
+ echo $i
+ newfile=`echo $i | sed 's/extra_nucl/nucl/g' `
+ echo $newfile
+ cp ${EZFIO_extra}/extra_nuclei/$i ${EZFIO_target}/nuclei/$newfile
+done
+
+### COPYING ALL DATA FROM $EZFIO_extra/ao_basis/ to ${EZFIO_target}/ao_extra_basis/
+direxists=false
+if [ -d ${EZFIO_target}/ao_extra_basis/ ] ; then
+ direxists=true
+ echo "The directory exists ao_extra_basis"
+ else
+ echo "PB !!"
+ direxists=false
+ mkdir ${EZFIO_target}/ao_extra_basis/
+fi
+echo "COPYING ALL DATA FROM "$EZFIO_extra"/ao_basis/ to "${EZFIO_target}"/ao_extra_basis/"
+data=`\ls $EZFIO_extra/ao_extra_basis/`
+for i in $data
+do
+ echo $i
+ newfile=`echo $i | sed 's/ao_extra/ao/g' `
+ echo $newfile
+ cp ${EZFIO_extra}/ao_extra_basis/$i ${EZFIO_target}/ao_basis/$newfile
+done
+
+
diff --git a/src/ao_extra_basis/qp_fit_1s_basis b/src/ao_extra_basis/qp_fit_1s_basis
new file mode 100755
index 00000000..366e478a
--- /dev/null
+++ b/src/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/src/ao_extra_basis/transform_basis.irp.f b/src/ao_extra_basis/transform_basis.irp.f
new file mode 100644
index 00000000..374d2dfc
--- /dev/null
+++ b/src/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/src/ao_extra_basis/uninstall b/src/ao_extra_basis/uninstall
new file mode 100755
index 00000000..0e4cd253
--- /dev/null
+++ b/src/ao_extra_basis/uninstall
@@ -0,0 +1,20 @@
+#!/bin/bash                                                                                                                                             
+
+# Check if the QP_ROOT environment variable is set. 
+if [[ -z ${QP_ROOT} ]]
+then
+  print "The QP_ROOT environment variable is not set."
+  print "Please reload the quantum_package.rc file."
+  exit -1
+fi
+
+# list of the scripts to be used by the module 
+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. 
+ for i in $scripts_list
+ do
+  find ${QP_ROOT}/scripts/$i -type l -delete
+ done
+
diff --git a/src/extra_nuclei/EZFIO.cfg b/src/extra_nuclei/EZFIO.cfg
new file mode 100644
index 00000000..8c04b4c7
--- /dev/null
+++ b/src/extra_nuclei/EZFIO.cfg
@@ -0,0 +1,56 @@
+[extra_nucl_num]
+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)
+size: (extra_nuclei.extra_nucl_num)
+interface: ezfio, provider
+
+[extra_nucl_charge]
+doc: Nuclear charges
+type:double precision
+size: (extra_nuclei.extra_nucl_num)
+interface: ezfio, provider
+
+[extra_nucl_coord]
+doc:  Nuclear coordinates in the format (:, {x,y,z})
+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/extra_nuclei/NEED b/src/extra_nuclei/NEED
new file mode 100644
index 00000000..0dc19a80
--- /dev/null
+++ b/src/extra_nuclei/NEED
@@ -0,0 +1,3 @@
+ezfio_files
+utils
+nuclei
diff --git a/src/extra_nuclei/README.rst b/src/extra_nuclei/README.rst
new file mode 100644
index 00000000..6bb260a0
--- /dev/null
+++ b/src/extra_nuclei/README.rst
@@ -0,0 +1,4 @@
+============
+extra_nuclei
+============
+
diff --git a/src/extra_nuclei/extra_nuclei.irp.f b/src/extra_nuclei/extra_nuclei.irp.f
new file mode 100644
index 00000000..53653e8d
--- /dev/null
+++ b/src/extra_nuclei/extra_nuclei.irp.f
@@ -0,0 +1,7 @@
+program extra_nuclei
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  print *, 'Hello world'
+end
diff --git a/src/extra_nuclei/nuclei.irp.f b/src/extra_nuclei/nuclei.irp.f
new file mode 100644
index 00000000..d8ea57f2
--- /dev/null
+++ b/src/extra_nuclei/nuclei.irp.f
@@ -0,0 +1,138 @@
+BEGIN_PROVIDER [ double precision, extra_nucl_coord,  (extra_nucl_num,3) ]
+   implicit none
+
+   BEGIN_DOC
+   ! Nuclear coordinates in the format (:, {x,y,z})
+   END_DOC
+   PROVIDE ezfio_filename extra_nucl_label extra_nucl_charge
+
+   if (mpi_master) then
+     double precision, allocatable  :: buffer(:,:)
+     extra_nucl_coord = 0.d0
+     allocate (buffer(extra_nucl_num,3))
+     buffer = 0.d0
+     logical                        :: has
+     call ezfio_has_extra_nuclei_extra_nucl_coord(has)
+     if (.not.has) then
+       print *, irp_here
+       stop 1
+     endif
+     call ezfio_get_extra_nuclei_extra_nucl_coord(buffer)
+     integer                        :: i,j
+
+     do i=1,3
+       do j=1,extra_nucl_num
+         extra_nucl_coord(j,i) = buffer(j,i)
+       enddo
+     enddo
+     deallocate(buffer)
+
+     character*(64), parameter      :: f = '(A16, 4(1X,F12.6))'
+     character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
+     double precision, parameter    :: a0= 0.529177249d0
+
+     call write_time(6)
+     write(6,'(A)') ''
+     write(6,'(A)') 'Extra Nuclear Coordinates (Angstroms)'
+     write(6,'(A)') '====================================='
+     write(6,'(A)') ''
+     write(6,ft)                                         &
+         '================','============','============','============','============'
+     write(6,*)                                          &
+         '     Atom          Charge          X            Y            Z '
+     write(6,ft)                                         &
+         '================','============','============','============','============'
+
+     do i=1,extra_nucl_num
+       write(6,f) extra_nucl_label(i), extra_nucl_charge(i),         &
+           extra_nucl_coord(i,1)*a0,                                       &
+           extra_nucl_coord(i,2)*a0,                                       &
+           extra_nucl_coord(i,3)*a0
+     enddo
+     write(6,ft)                                         &
+         '================','============','============','============','============'
+     write(6,'(A)') ''
+
+     if (extra_nucl_num > 1) then
+       double precision               :: dist_min, x, y, z
+       dist_min = huge(1.d0)
+       do i=1,extra_nucl_num
+         do j=i+1,extra_nucl_num
+           x = extra_nucl_coord(i,1)-extra_nucl_coord(j,1)
+           y = extra_nucl_coord(i,2)-extra_nucl_coord(j,2)
+           z = extra_nucl_coord(i,3)-extra_nucl_coord(j,3)
+           dist_min = min(x*x + y*y + z*z, dist_min)
+         enddo
+       enddo
+       write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', &
+          dsqrt(dist_min)*a0,' Angstrom'
+     endif
+
+   endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+   IRP_IF MPI
+     include 'mpif.h'
+     integer                        :: ierr
+     call MPI_BCAST( extra_nucl_coord, 3*extra_nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+     if (ierr /= MPI_SUCCESS) then
+       stop 'Unable to read nucl_coord with MPI'
+     endif
+   IRP_ENDIF
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, extra_nucl_coord_transp, (3,extra_nucl_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Transposed array of extra_nucl_coord
+   END_DOC
+   integer                        :: i, k
+   extra_nucl_coord_transp = 0.d0
+
+   do i=1,extra_nucl_num
+     extra_nucl_coord_transp(1,i) = extra_nucl_coord(i,1)
+     extra_nucl_coord_transp(2,i) = extra_nucl_coord(i,2)
+     extra_nucl_coord_transp(3,i) = extra_nucl_coord(i,3)
+   enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, extra_center_of_mass, (3) ]
+  implicit none
+  BEGIN_DOC
+  ! Center of mass of the molecule
+  END_DOC
+  integer                        :: i,j
+  double precision               :: s
+  extra_center_of_mass(:) = 0.d0
+  s = 0.d0
+  do i=1,extra_nucl_num
+    do j=1,3
+      extra_center_of_mass(j) += extra_nucl_coord(i,j)* element_mass(int(extra_nucl_charge(i)))
+    enddo
+    s += element_mass(int(extra_nucl_charge(i)))
+  enddo
+  s = 1.d0/s
+  extra_center_of_mass(:) = extra_center_of_mass(:)*s
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, extra_nucl_dist, (extra_nucl_num,extra_nucl_num)]
+ implicit none
+ integer :: i,j
+ double precision :: x,y,z
+ do i = 1, extra_nucl_num
+  do j = 1, extra_nucl_num
+   x = extra_nucl_coord(i,1)-extra_nucl_coord(j,1)
+   y = extra_nucl_coord(i,2)-extra_nucl_coord(j,2)
+   z = extra_nucl_coord(i,3)-extra_nucl_coord(j,3)
+   extra_nucl_dist(j,i) = x*x+y*y+z*z
+   extra_nucl_dist(j,i) = dsqrt(extra_nucl_dist(j,i))
+  enddo
+ enddo
+
+END_PROVIDER 
diff --git a/src/mo_one_e_ints/mo_one_e_ints.irp.f b/src/mo_one_e_ints/mo_one_e_ints.irp.f
index 19cd2159..db05ba90 100644
--- a/src/mo_one_e_ints/mo_one_e_ints.irp.f
+++ b/src/mo_one_e_ints/mo_one_e_ints.irp.f
@@ -22,13 +22,3 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)]
 
 END_PROVIDER
 
-
-BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)]
- implicit none
- BEGIN_DOC
-! Integrals of the one e hamiltonian obtained from the integrals on the MO basis 
-!
-! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions
- END_DOC
- call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num)
-END_PROVIDER