diff --git a/external/irpf90 b/external/irpf90
index 43160c60..4ab1b175 160000
--- a/external/irpf90
+++ b/external/irpf90
@@ -1 +1 @@
-Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b
+Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
diff --git a/plugins/local/extra_basis_int/ao_overlap.irp.f b/plugins/local/extra_basis_int/ao_overlap.irp.f
index 4f8debb6..9e45e690 100644
--- a/plugins/local/extra_basis_int/ao_overlap.irp.f
+++ b/plugins/local/extra_basis_int/ao_overlap.irp.f
@@ -4,7 +4,7 @@
  BEGIN_PROVIDER [double precision, ao_extra_overlap  , (ao_extra_num, ao_extra_num)]
 
   BEGIN_DOC
-  ! Overlap between atomic basis functions:
+  ! Overlap between atomic basis functions belonging to the EXTRA BASIS 
   !
   ! :math:`\int \chi_i(r) \chi_j(r) dr`
   END_DOC
@@ -69,7 +69,9 @@ END_PROVIDER
  BEGIN_PROVIDER [double precision, ao_extra_overlap_mixed  , (ao_num, ao_extra_num)]
 
   BEGIN_DOC
-  ! Overlap between atomic basis functions:
+  ! Overlap between atomic basis functions:  
+  !
+  ! first index belongs to the REGULAR AO basis, second to the EXTRA basis 
   !
   ! <AO_i|AO_j extra basis>
   END_DOC
diff --git a/plugins/local/extra_basis_int/coul_1s.irp.f b/plugins/local/extra_basis_int/coul_1s.irp.f
index 964ed36a..62ec4331 100644
--- a/plugins/local/extra_basis_int/coul_1s.irp.f
+++ b/plugins/local/extra_basis_int/coul_1s.irp.f
@@ -11,6 +11,11 @@ double precision function coul_full_ao_pq_r_1s(p,q,R,R_p,R_q)
  double precision, intent(in) :: R(3),R_p(3),R_q(3)
  integer, intent(in) :: p,q
  double precision :: coef,dist,P_pq(3),coefaos
+ if(.not.ao_extra_only_1s)then
+   print*,'You are using a function assuming that the extra basis is fitted on 1s functions'
+   print*,'But this is not the case apparently ... stopping'
+   stop
+ endif
  coefaos= ao_extra_coef_normalized(p,1) * ao_extra_coef_normalized(q,1) 
  coef = inv_pi_gamma_pq_3_2_ao_extra(p,q) * E_pq_ao_extra(p,q)
  P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q
@@ -40,6 +45,11 @@ double precision function coul_pq_r_1s(p,q,R,R_p,R_q)
  double precision, intent(in) :: R(3),R_p(3),R_q(3)
  integer, intent(in) :: p,q
  double precision :: dist,P_pq(3)
+ if(.not.ao_extra_only_1s)then
+   print*,'You are using a function assuming that the extra basis is fitted on 1s functions'
+   print*,'But this is not the case apparently ... stopping'
+   stop
+ endif
  P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q
  P_pq = P_pq * inv_gamma_pq_ao_extra(q,p)
  dist = (P_pq(1)-R(1)) * (P_pq(1)-R(1))
diff --git a/plugins/local/extra_basis_int/extra_basis_int.irp.f b/plugins/local/extra_basis_int/extra_basis_int.irp.f
index 1d35b1c2..59b90374 100644
--- a/plugins/local/extra_basis_int/extra_basis_int.irp.f
+++ b/plugins/local/extra_basis_int/extra_basis_int.irp.f
@@ -11,18 +11,30 @@ program extra_basis_int
 ! call routine_pot_ne
 ! call routine_test_pot_ne_extra_mixed
 ! call routine_test_coul_1s
- call print_v_ne_extra_basis
- call print_v_ne_basis
+! call print_v_ne_extra_basis
+! call print_v_ne_basis
+! call test_v_ne_a_extra_basis
+! call print_v_ee_mixed_direct
+ call print_v_ee_mixed_exchange
  
 end
 
+subroutine test_v_ne_a_extra_basis
+ implicit none
+ integer :: i,j
+ do i = 1, ao_extra_num
+  write(*,'(100(F16.10,X))')pot_vne_A_extra_basis(1:ao_extra_num,i)
+ enddo
+end
+
+
 subroutine test_overlap
  implicit none
   integer :: i,j
-  do i = 1, ao_extra_num
-   do j = 1, ao_extra_num 
-    write(33,*)ao_extra_overlap(j,i)
-   enddo
+  do i = 1, ao_num
+!   do j = 1, ao_num 
+    write(33,'(100(F16.10,X))')ao_extra_overlap_mixed(i,1:ao_extra_num)
+!   enddo
   enddo
 end
 
@@ -189,3 +201,35 @@ subroutine print_v_ne_basis
  print*,'accu = ',accu
 
 end
+
+subroutine print_v_ee_mixed_direct
+ implicit none
+ integer :: i,j,k,l
+ double precision :: ao_two_e_integral_mixed_direct
+ do i = 1, ao_num
+  do j = 1, ao_num
+   do k = 1, ao_extra_num
+    do l = 1, ao_extra_num
+     write(34,*)ao_two_e_integral_mixed_direct(i, j, k, l)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine print_v_ee_mixed_exchange
+ implicit none
+ integer :: i,j,k,l
+ double precision :: ao_two_e_integral_mixed_exchange
+ do i = 1, ao_num
+  do j = 1, ao_extra_num
+   do k = 1, ao_num
+    do l = 1, ao_extra_num
+     write(34,*)ao_two_e_integral_mixed_exchange(i, j, k, l)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
diff --git a/plugins/local/extra_basis_int/pot_ao_ints.irp.f b/plugins/local/extra_basis_int/pot_ao_ints.irp.f
index 5f3af244..f52e7579 100644
--- a/plugins/local/extra_basis_int/pot_ao_ints.irp.f
+++ b/plugins/local/extra_basis_int/pot_ao_ints.irp.f
@@ -1,3 +1,53 @@
+BEGIN_PROVIDER [ double precision, pot_vne_A_extra_basis, (ao_extra_num,ao_extra_num)]
+ implicit none
+  BEGIN_DOC
+  !
+  ! Computes the following integral :
+  ! $\sum_{R in the USUAL nuclei} -Z <chi_i|1/|r-R||chi_j>$ 
+  !
+  ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis 
+  END_DOC
+  integer :: mu,nu
+  double precision :: v_nucl_extra_ao
+  pot_vne_A_extra_basis = 0.d0
+  do mu = 1, ao_extra_num
+   do nu = 1, ao_extra_num
+    pot_vne_A_extra_basis(nu,mu)= v_nucl_extra_ao(mu,nu)
+   enddo
+  enddo
+
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, pot_vne_extra_basis, (ao_num,ao_num)]
+ implicit none
+  BEGIN_DOC
+  !
+  ! Computes the following integral :
+  ! $\sum_{R in EXTRA nuclei} -Z <chi_i|1/|r-R||chi_j>$ 
+  !
+  !
+  ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the USUAL basis 
+  END_DOC
+  integer :: mu,nu,k_nucl
+  double precision :: mu_in, R_nucl(3),charge_nucl, integral
+  double precision :: NAI_pol_mult_erf_ao
+  mu_in = 10.d0**10
+  pot_vne_extra_basis = 0.d0
+  do mu = 1, ao_num
+   do nu = 1, ao_num
+    do k_nucl = 1, extra_nucl_num
+     R_nucl(1:3) = extra_nucl_coord_transp(1:3,k_nucl)
+     charge_nucl = extra_nucl_charge(k_nucl)
+     integral = NAI_pol_mult_erf_ao(mu, nu, mu_in, R_nucl)
+     pot_vne_extra_basis(nu,mu) += -integral * charge_nucl
+    enddo
+   enddo
+  enddo
+
+END_PROVIDER 
+
+
 
 double precision function NAI_pol_mult_erf_ao_extra(i_ao, j_ao, mu_in, C_center)
 
diff --git a/plugins/local/extra_basis_int/ref_extra_basis.irp.f b/plugins/local/extra_basis_int/ref_extra_basis.irp.f
index 39055fd0..70d77733 100644
--- a/plugins/local/extra_basis_int/ref_extra_basis.irp.f
+++ b/plugins/local/extra_basis_int/ref_extra_basis.irp.f
@@ -6,7 +6,9 @@ program pouet
 ! call routine_pot_ne_extra
 ! call ref_pot_ne_mixed
 ! call ref_pot_ne
- call ref_pot_ne_extra_mixed
+! call ref_pot_ne_extra_mixed
+! call ref_v_ee_mixed_direct
+ call ref_v_ee_mixed_exchange
 
 end
 
@@ -113,3 +115,35 @@ subroutine ref_pot_ne_extra_mixed
   enddo
  enddo
 end
+
+subroutine ref_v_ee_mixed_direct
+ implicit none
+ integer :: i,j,k,l
+ double precision :: ao_two_e_integral
+ do i = 1, 15
+  do j = 1, 15
+   do k = 16, ao_num
+    do l = 16, ao_num
+     write(33,*)ao_two_e_integral(i, j, k, l)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine ref_v_ee_mixed_exchange
+ implicit none
+ integer :: i,j,k,l
+ double precision :: ao_two_e_integral
+ do i = 1, 15
+  do j = 16, ao_num
+   do k = 1, 15
+    do l = 16, ao_num
+     write(33,*)ao_two_e_integral(i, j, k, l)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
diff --git a/plugins/local/extra_basis_int/two_e_int.irp.f b/plugins/local/extra_basis_int/two_e_int.irp.f
new file mode 100644
index 00000000..2cde4153
--- /dev/null
+++ b/plugins/local/extra_basis_int/two_e_int.irp.f
@@ -0,0 +1,145 @@
+double precision function ao_two_e_integral_mixed_direct(i, j, k, l)
+
+  BEGIN_DOC
+  !  integral of the AO basis <ik|jl> or (ij|kl)
+  !     i(r1) j(r1) 1/r12 k(r2) l(r2)
+  !     A     A           B     B
+  !
+  ! where i,j belong to the REGULAR AO basis (system A) and k,l to the EXTRA basis (system B)
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer, intent(in)            :: i, j, k, l
+
+  integer                        :: p, q, r, s
+  integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
+  integer                        :: iorder_p(3), iorder_q(3)
+  double precision               :: I_center(3), J_center(3), K_center(3), L_center(3)
+  double precision               :: integral
+  double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
+  double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
+  double precision               :: general_primitive_integral
+
+  dim1 = n_pt_max_integrals
+
+  num_i = ao_nucl(i)
+  num_j = ao_nucl(j)
+  num_k = ao_extra_nucl(k)
+  num_l = ao_extra_nucl(l)
+  ao_two_e_integral_mixed_direct = 0.d0
+
+  do p = 1, 3
+    I_power(p) = ao_power(i,p)
+    J_power(p) = ao_power(j,p)
+    K_power(p) = ao_extra_power(k,p)
+    L_power(p) = ao_extra_power(l,p)
+    I_center(p) = nucl_coord(num_i,p)
+    J_center(p) = nucl_coord(num_j,p)
+    K_center(p) = extra_nucl_coord(num_k,p)
+    L_center(p) = extra_nucl_coord(num_l,p)
+  enddo
+
+  double precision               :: coef1, coef2, coef3, coef4
+  double precision               :: p_inv,q_inv
+
+  do p = 1, ao_prim_num(i)
+    coef1 = ao_coef_normalized_ordered_transp(p,i)
+    do q = 1, ao_prim_num(j)
+      coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
+      call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
+          ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),                 &
+          I_power,J_power,I_center,J_center,dim1)
+      p_inv = 1.d0/pp
+      do r = 1, ao_extra_prim_num(k)
+        coef3 = coef2*ao_extra_coef_normalized_ordered_transp(r,k)
+        do s = 1, ao_extra_prim_num(l)
+          coef4 = coef3*ao_extra_coef_normalized_ordered_transp(s,l)
+          call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
+              ao_extra_expo_ordered_transp(r,k),ao_extra_expo_ordered_transp(s,l),             &
+              K_power,L_power,K_center,L_center,dim1)
+          q_inv = 1.d0/qq
+          integral = general_primitive_integral(dim1,              &
+              P_new,P_center,fact_p,pp,p_inv,iorder_p,             &
+              Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
+          ao_two_e_integral_mixed_direct = ao_two_e_integral_mixed_direct +  coef4 * integral
+        enddo ! s
+      enddo  ! r
+    enddo   ! q
+  enddo    ! p
+
+end
+
+double precision function ao_two_e_integral_mixed_exchange(i, j, k, l)
+
+  BEGIN_DOC
+  !  integral of the AO basis <ik|jl> or (ij|kl)
+  !     i(r1) j(r1) 1/r12 k(r2) l(r2)
+  !     A     B           A     B
+  !
+  ! where i,k belong to the REGULAR AO basis (system A) and j,l to the EXTRA basis (system B)
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer, intent(in)            :: i, j, k, l
+
+  integer                        :: p, q, r, s
+  integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
+  integer                        :: iorder_p(3), iorder_q(3)
+  double precision               :: I_center(3), J_center(3), K_center(3), L_center(3)
+  double precision               :: integral
+  double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
+  double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
+  double precision               :: general_primitive_integral
+
+  dim1 = n_pt_max_integrals
+
+  num_i = ao_nucl(i)
+  num_j = ao_extra_nucl(j)
+  num_k = ao_nucl(k)
+  num_l = ao_extra_nucl(l)
+  ao_two_e_integral_mixed_exchange = 0.d0
+
+  do p = 1, 3
+    I_power(p) = ao_power(i,p)
+    J_power(p) = ao_extra_power(j,p)
+    K_power(p) = ao_power(k,p)
+    L_power(p) = ao_extra_power(l,p)
+    I_center(p) = nucl_coord(num_i,p)
+    J_center(p) = extra_nucl_coord(num_j,p)
+    K_center(p) = nucl_coord(num_k,p)
+    L_center(p) = extra_nucl_coord(num_l,p)
+  enddo
+
+  double precision               :: coef1, coef2, coef3, coef4
+  double precision               :: p_inv,q_inv
+
+  do p = 1, ao_prim_num(i)
+    coef1 = ao_coef_normalized_ordered_transp(p,i)
+    do q = 1, ao_extra_prim_num(j)
+      coef2 = coef1*ao_extra_coef_normalized_ordered_transp(q,j)
+      call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
+          ao_expo_ordered_transp(p,i),ao_extra_expo_ordered_transp(q,j),                 &
+          I_power,J_power,I_center,J_center,dim1)
+      p_inv = 1.d0/pp
+      do r = 1, ao_prim_num(k)
+        coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
+        do s = 1, ao_extra_prim_num(l)
+          coef4 = coef3*ao_extra_coef_normalized_ordered_transp(s,l)
+          call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
+              ao_expo_ordered_transp(r,k),ao_extra_expo_ordered_transp(s,l),             &
+              K_power,L_power,K_center,L_center,dim1)
+          q_inv = 1.d0/qq
+          integral = general_primitive_integral(dim1,              &
+              P_new,P_center,fact_p,pp,p_inv,iorder_p,             &
+              Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
+          ao_two_e_integral_mixed_exchange = ao_two_e_integral_mixed_exchange +  coef4 * integral
+        enddo ! s
+      enddo  ! r
+    enddo   ! q
+  enddo    ! p
+
+end
diff --git a/plugins/local/extra_basis_int/v_mixed_extra.irp.f b/plugins/local/extra_basis_int/v_mixed_extra.irp.f
index ac856692..8b8ce92e 100644
--- a/plugins/local/extra_basis_int/v_mixed_extra.irp.f
+++ b/plugins/local/extra_basis_int/v_mixed_extra.irp.f
@@ -1,3 +1,6 @@
+!!! TODO:: optimize when "ao_extra_only_1s" is True
+
+
 double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
  implicit none
   BEGIN_DOC
@@ -6,9 +9,9 @@ double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
   ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne^{extra}(r)$.
   !
   !
-  ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis 
+  ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis (system B)
   !
-  ! and v_ne^{extra}(r) is the Coulomb potential coming from the EXTRA nuclei
+  ! and v_ne^{extra}(r) is the Coulomb potential coming from the EXTRA nuclei (system B)
   END_DOC
  integer, intent(in) ::i_ao,j_ao
  double precision :: mu_in,charge,coord(3)
@@ -23,6 +26,30 @@ double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
  enddo
 end
 
+double precision function v_extra_nucl_ao(i_ao,j_ao)
+ implicit none
+  BEGIN_DOC
+  !
+  ! Computes the following integral :
+  ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne(r)$.
+  !
+  !
+  ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the REGULAR basis (system A)
+  !
+  ! and v_ne(r) is the Coulomb potential coming from the EXTRA nuclei (system B)
+  END_DOC
+ integer, intent(in) ::i_ao,j_ao
+ integer :: i
+ double precision :: mu_in, coord(3),charge, integral
+ double precision :: NAI_pol_mult_erf_ao
+ mu_in = 1.d+10
+ do i = 1, extra_nucl_num
+  coord(1:3) = extra_nucl_coord_transp(1:3,i)
+  charge = extra_nucl_charge(i)
+  v_extra_nucl_ao += -NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, coord) * charge
+ enddo
+end
+
 
 double precision function v_nucl_extra_ao(i_ao,j_ao)
  implicit none
@@ -32,9 +59,9 @@ double precision function v_nucl_extra_ao(i_ao,j_ao)
   ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne(r)$.
   !
   !
-  ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis 
+  ! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis (system B)
   !
-  ! and v_ne(r) is the Coulomb potential coming from the REGULAR nuclei
+  ! and v_ne(r) is the Coulomb potential coming from the REGULAR nuclei (system A)
   END_DOC
  integer, intent(in) ::i_ao,j_ao
  double precision :: mu_in,charge,coord(3)
diff --git a/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f b/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f
index e71d49fc..68d965cb 100644
--- a/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f
+++ b/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f
@@ -30,6 +30,7 @@ BEGIN_PROVIDER [ double precision, trace_ao_one_e_ints]
 !          have the same number of functions
  END_DOC
  integer :: i,j
+ double precision :: accu
  double precision, allocatable :: inv_overlap_times_integrals(:,:) ! = h S^{-1}
  allocate(inv_overlap_times_integrals(ao_num,ao_num))
  ! routine that computes the product of two matrices, you can check it with 
diff --git a/src/ao_extra_basis/README.rst b/src/ao_extra_basis/README.rst
index f60d71c0..5f850255 100644
--- a/src/ao_extra_basis/README.rst
+++ b/src/ao_extra_basis/README.rst
@@ -5,7 +5,9 @@ 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. 
+Check in the directory "tuto" for a simple example of how to create a fictious system "B" attached independently to a system "A"
+
+Another 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
@@ -13,3 +15,4 @@ This is done with the various scripts here:
  
 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/fit_1s_basis.irp.f b/src/ao_extra_basis/fit_1s_basis.irp.f
index ef09d5b2..9c788173 100644
--- a/src/ao_extra_basis/fit_1s_basis.irp.f
+++ b/src/ao_extra_basis/fit_1s_basis.irp.f
@@ -31,6 +31,7 @@ program fit_1s_basis
  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_only_1s(.True.)
  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)
diff --git a/src/ao_extra_basis/qp_copy_extra_basis b/src/ao_extra_basis/qp_copy_extra_basis
index cb435e18..6d0e17b8 100755
--- a/src/ao_extra_basis/qp_copy_extra_basis
+++ b/src/ao_extra_basis/qp_copy_extra_basis
@@ -58,7 +58,7 @@ do
 done
 i=primitives_normalized
 newfile=primitives_normalized_extra
-cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
+cp ${EZFIO_extra}/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
diff --git a/src/ao_extra_basis/tuto/He_A.xyz b/src/ao_extra_basis/tuto/He_A.xyz
new file mode 100644
index 00000000..d5285ade
--- /dev/null
+++ b/src/ao_extra_basis/tuto/He_A.xyz
@@ -0,0 +1,3 @@
+1
+He atom "A"
+He 0. 0. 0.
diff --git a/src/ao_extra_basis/tuto/example_copy.sh b/src/ao_extra_basis/tuto/example_copy.sh
new file mode 100755
index 00000000..0677b183
--- /dev/null
+++ b/src/ao_extra_basis/tuto/example_copy.sh
@@ -0,0 +1,26 @@
+source ~/qp2/quantum_package.rc
+## Example of how to generate an additional h2o molecule, stored as a extra basis/nuclei etc .. to an He
+
+sys_B=h2o.xyz
+basis_B=sto-3g
+output_B=${sys_B%.xyz}_${basis_B}
+
+sys_A=He_A.xyz
+basis_A=cc-pvtz
+output_A=${sys_A%.xyz}_${basis_A}_extra_${output_B}
+
+# we create the system "B" that will be attached as an "extra system" to the syste "A"
+qp create_ezfio -b $basis_B $sys_B -o ${output_B}
+# we perform an HF calculation to obtain the AO density matrix 
+qp run scf 
+# we save the density matrix in the EZFIO
+qp run save_one_e_dm 
+# we create the system "A" 
+qp create_ezfio -b $basis_A $sys_A -o ${output_A}
+# We perform an SCF calculation
+qp run scf 
+# we copy the system "B" information as extra nuclei/basis etc in the EZFIO of system "A"
+qp_copy_extra_basis ${output_B} ${output_A}
+
+# we execute an example of progra that prints a lot of useful integrals/information on the A-B interaction
+qp run test_extra_basis | tee ${output_A}.test_extra_basis
diff --git a/src/ao_extra_basis/tuto/h2o.xyz b/src/ao_extra_basis/tuto/h2o.xyz
new file mode 100644
index 00000000..d3928214
--- /dev/null
+++ b/src/ao_extra_basis/tuto/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   
+
+