diff --git a/plugins/local/ao_extra_basis/EZFIO.cfg b/plugins/local/ao_extra_basis/EZFIO.cfg
deleted file mode 100644
index 79b37bf0..00000000
--- a/plugins/local/ao_extra_basis/EZFIO.cfg
+++ /dev/null
@@ -1,104 +0,0 @@
-[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,1)
-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
deleted file mode 100644
index 03431803..00000000
--- a/plugins/local/ao_extra_basis/LiH.xyz
+++ /dev/null
@@ -1,4 +0,0 @@
-2
-
-H 0. 0. 0.
-Li 0. 0. 1.0
diff --git a/plugins/local/ao_extra_basis/NEED b/plugins/local/ao_extra_basis/NEED
deleted file mode 100644
index c4e8c3cf..00000000
--- a/plugins/local/ao_extra_basis/NEED
+++ /dev/null
@@ -1,3 +0,0 @@
-extra_nuclei
-basis
-ao_basis
diff --git a/plugins/local/ao_extra_basis/README.rst b/plugins/local/ao_extra_basis/README.rst
deleted file mode 100644
index f60d71c0..00000000
--- a/plugins/local/ao_extra_basis/README.rst
+++ /dev/null
@@ -1,15 +0,0 @@
-===========
-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/aos.irp.f b/plugins/local/ao_extra_basis/aos.irp.f
deleted file mode 100644
index 56d6fb04..00000000
--- a/plugins/local/ao_extra_basis/aos.irp.f
+++ /dev/null
@@ -1,325 +0,0 @@
-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/plugins/local/ao_extra_basis/aos_transp.irp.f b/plugins/local/ao_extra_basis/aos_transp.irp.f
deleted file mode 100644
index ed34835b..00000000
--- a/plugins/local/ao_extra_basis/aos_transp.irp.f
+++ /dev/null
@@ -1,68 +0,0 @@
-
-! ---
-
-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/plugins/local/ao_extra_basis/density_extra.irp.f b/plugins/local/ao_extra_basis/density_extra.irp.f
deleted file mode 100644
index fa215466..00000000
--- a/plugins/local/ao_extra_basis/density_extra.irp.f
+++ /dev/null
@@ -1,87 +0,0 @@
-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,1) * 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/plugins/local/ao_extra_basis/dimensions_integrals.irp.f b/plugins/local/ao_extra_basis/dimensions_integrals.irp.f
deleted file mode 100644
index 80dcda18..00000000
--- a/plugins/local/ao_extra_basis/dimensions_integrals.irp.f
+++ /dev/null
@@ -1,19 +0,0 @@
- 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/plugins/local/ao_extra_basis/extra_basis.irp.f b/plugins/local/ao_extra_basis/extra_basis.irp.f
deleted file mode 100644
index 9773af2f..00000000
--- a/plugins/local/ao_extra_basis/extra_basis.irp.f
+++ /dev/null
@@ -1,16 +0,0 @@
-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/plugins/local/ao_extra_basis/fit_1s_basis.irp.f b/plugins/local/ao_extra_basis/fit_1s_basis.irp.f
deleted file mode 100644
index ef09d5b2..00000000
--- a/plugins/local/ao_extra_basis/fit_1s_basis.irp.f
+++ /dev/null
@@ -1,41 +0,0 @@
-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/plugins/local/ao_extra_basis/h2o.xyz b/plugins/local/ao_extra_basis/h2o.xyz
deleted file mode 100644
index d3928214..00000000
--- a/plugins/local/ao_extra_basis/h2o.xyz
+++ /dev/null
@@ -1,7 +0,0 @@
-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
deleted file mode 100755
index e7e668e1..00000000
--- a/plugins/local/ao_extra_basis/install
+++ /dev/null
@@ -1,23 +0,0 @@
-#!/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/plugins/local/ao_extra_basis/prov_fit_1s.irp.f b/plugins/local/ao_extra_basis/prov_fit_1s.irp.f
deleted file mode 100644
index 99a18a6b..00000000
--- a/plugins/local/ao_extra_basis/prov_fit_1s.irp.f
+++ /dev/null
@@ -1,266 +0,0 @@
-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/plugins/local/ao_extra_basis/qp_add_extra_fit_system b/plugins/local/ao_extra_basis/qp_add_extra_fit_system
deleted file mode 100755
index dae95062..00000000
--- a/plugins/local/ao_extra_basis/qp_add_extra_fit_system
+++ /dev/null
@@ -1,19 +0,0 @@
-#!/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
deleted file mode 100755
index cb435e18..00000000
--- a/plugins/local/ao_extra_basis/qp_copy_extra_basis
+++ /dev/null
@@ -1,67 +0,0 @@
-#!/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/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
deleted file mode 100755
index ba7e6a71..00000000
--- a/plugins/local/ao_extra_basis/qp_copy_extra_basis_to_usual_basis
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/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/plugins/local/ao_extra_basis/qp_fit_1s_basis b/plugins/local/ao_extra_basis/qp_fit_1s_basis
deleted file mode 100755
index 366e478a..00000000
--- a/plugins/local/ao_extra_basis/qp_fit_1s_basis
+++ /dev/null
@@ -1,12 +0,0 @@
-#!/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
deleted file mode 100644
index 374d2dfc..00000000
--- a/plugins/local/ao_extra_basis/transform_basis.irp.f
+++ /dev/null
@@ -1,49 +0,0 @@
-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
deleted file mode 100755
index 0e4cd253..00000000
--- a/plugins/local/ao_extra_basis/uninstall
+++ /dev/null
@@ -1,20 +0,0 @@
-#!/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/plugins/local/extra_nuclei/EZFIO.cfg b/plugins/local/extra_nuclei/EZFIO.cfg
deleted file mode 100644
index 8c04b4c7..00000000
--- a/plugins/local/extra_nuclei/EZFIO.cfg
+++ /dev/null
@@ -1,56 +0,0 @@
-[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/plugins/local/extra_nuclei/NEED b/plugins/local/extra_nuclei/NEED
deleted file mode 100644
index 0dc19a80..00000000
--- a/plugins/local/extra_nuclei/NEED
+++ /dev/null
@@ -1,3 +0,0 @@
-ezfio_files
-utils
-nuclei
diff --git a/plugins/local/extra_nuclei/README.rst b/plugins/local/extra_nuclei/README.rst
deleted file mode 100644
index 6bb260a0..00000000
--- a/plugins/local/extra_nuclei/README.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-============
-extra_nuclei
-============
-
diff --git a/plugins/local/extra_nuclei/extra_nuclei.irp.f b/plugins/local/extra_nuclei/extra_nuclei.irp.f
deleted file mode 100644
index 53653e8d..00000000
--- a/plugins/local/extra_nuclei/extra_nuclei.irp.f
+++ /dev/null
@@ -1,7 +0,0 @@
-program extra_nuclei
-  implicit none
-  BEGIN_DOC
-! TODO : Put the documentation of the program here
-  END_DOC
-  print *, 'Hello world'
-end
diff --git a/plugins/local/extra_nuclei/nuclei.irp.f b/plugins/local/extra_nuclei/nuclei.irp.f
deleted file mode 100644
index d8ea57f2..00000000
--- a/plugins/local/extra_nuclei/nuclei.irp.f
+++ /dev/null
@@ -1,138 +0,0 @@
-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/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
index 69fa16ff..6e9d58fd 100644
--- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f
+++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
@@ -122,7 +122,7 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
 
  ! Contracted two-e integrals : intermediate quantity
  ! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n >
- allocate( v_tilde(n_act_orb,n_act_orb)   )
+ allocate( v_tilde(n_inact_orb,n_act_orb)   )
  allocate( integrals_array(mo_num,mo_num) )
  v_tilde = 0.d0
  do a = 1, n_act_orb