diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py
index a515efba..b8d0d8ed 100755
--- a/scripts/qp_import_trexio.py
+++ b/scripts/qp_import_trexio.py
@@ -271,12 +271,11 @@ def write_ezfio(trexio_filename, filename):
     if basis_type.lower() == "gaussian" and not cartesian:
         try:
           import trexio_tools
-          fd, tmp = tempfile.mkstemp()
-          os.close(fd)
+          tmp = "cartesian_"+trexio_filename
           retcode = subprocess.call(["trexio", "convert-to", "-t", "cartesian", "-o", tmp, trexio_filename])
           trexio_file_cart = trexio.File(tmp,mode='r',back_end=trexio.TREXIO_AUTO)
           cartesian = trexio.read_ao_cartesian(trexio_file_cart)
-          os.unlink(tmp)
+          ezfio.set_trexio_trexio_file(tmp)
         except:
           pass
 
diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f
index 440cc865..d9381015 100644
--- a/src/ao_basis/aos.irp.f
+++ b/src/ao_basis/aos.irp.f
@@ -20,7 +20,22 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ]
      ao_shell(k) = i
    enddo
  enddo
+END_PROVIDER
 
+BEGIN_PROVIDER [ integer, ao_sphe_shell, (ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the AO corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)
+   do j=-n,n
+     k = k+1
+     ao_sphe_shell(k) = i
+   enddo
+ enddo
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer, ao_first_of_shell, (shell_num) ]
@@ -133,6 +148,16 @@ END_PROVIDER
     enddo
   enddo
 
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_sphe_coef_normalization_factor, (ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Normalization factor in spherical AO basis
+ END_DOC
+
+ ao_sphe_coef_normalization_factor(:) = 1.d0
+
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num,ao_prim_num_max) ]
diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f
index 336161f8..53a60413 100644
--- a/src/ao_basis/spherical_to_cartesian.irp.f
+++ b/src/ao_basis/spherical_to_cartesian.irp.f
@@ -4,7 +4,8 @@
 ! First index is the index of the cartesian AO, obtained by ao_power_index
 ! Second index is the index of the spherical AO
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_0, (1) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=0
@@ -12,10 +13,12 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ]
   cart_to_sphe_0 = 0.d0
 
   cart_to_sphe_0 ( 1, 1) = 1.0d0
+  cart_to_sphe_norm_0 (1) = 1.d0
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_1, (3) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=1
@@ -25,10 +28,14 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ]
   cart_to_sphe_1 ( 3, 1) = 1.0d0
   cart_to_sphe_1 ( 1, 2) = 1.0d0
   cart_to_sphe_1 ( 2, 3) = 1.0d0
+  cart_to_sphe_norm_1 (1) = 1.d0
+  cart_to_sphe_norm_1 (2) = 1.d0
+  cart_to_sphe_norm_1 (3) = 1.d0
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_2, (6) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=2
@@ -43,10 +50,14 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ]
   cart_to_sphe_2 ( 1, 4) = 0.86602540378443864676d0
   cart_to_sphe_2 ( 4, 4) = -0.86602540378443864676d0
   cart_to_sphe_2 ( 2, 5) = 1.0d0
+
+  cart_to_sphe_norm_2 = (/ 1.0d0, 1.7320508075688772d0, 1.7320508075688772d0, 1.0d0, &
+                           1.7320508075688772d0, 1.0d0 /)
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_3, (10) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=3
@@ -69,10 +80,15 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ]
   cart_to_sphe_3 ( 4, 6) = -1.0606601717798212866d0
   cart_to_sphe_3 ( 2, 7) = 1.0606601717798212866d0
   cart_to_sphe_3 ( 7, 7) = -0.790569415042094833d0
+
+  cart_to_sphe_norm_3 = (/ 1.0d0, 2.23606797749979d0, 2.23606797749979d0, &
+  2.23606797749979d0, 3.872983346207417d0, 2.23606797749979d0, 1.0d0, 2.23606797749979d0, &
+  2.23606797749979d0, 1.d00 /)
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_4, (15) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=4
@@ -107,10 +123,18 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ]
   cart_to_sphe_4 (11, 8) = 0.73950997288745200532d0
   cart_to_sphe_4 ( 2, 9) = 1.1180339887498948482d0
   cart_to_sphe_4 ( 7, 9) = -1.1180339887498948482d0
+
+  cart_to_sphe_norm_4 = (/ 1.0d0, 2.6457513110645907d0, 2.6457513110645907d0, &
+  3.4156502553198664d0, 5.916079783099616d0, 3.415650255319866d0, &
+  2.6457513110645907d0, 5.916079783099616d0, 5.916079783099616d0, &
+  2.6457513110645907d0, 1.0d0, 2.6457513110645907d0, 3.415650255319866d0, &
+  2.6457513110645907d0, 1.d00 /)
+
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_5, (21) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=5
@@ -163,10 +187,18 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ]
   cart_to_sphe_5 ( 2,11) = 1.169267933366856683d0
   cart_to_sphe_5 ( 7,11) = -1.5309310892394863114d0
   cart_to_sphe_5 (16,11) = 0.7015607600201140098d0
+
+  cart_to_sphe_norm_5 = (/ 1.0d0, 3.0d0, 3.0d0, 4.58257569495584d0, &
+  7.937253933193773d0, 4.58257569495584d0, 4.58257569495584d0, &
+  10.246950765959598d0, 10.246950765959598d0, 4.582575694955841d0, 3.0d0, &
+  7.937253933193773d0, 10.246950765959598d0, 7.937253933193773d0, 3.0d0, 1.0d0, &
+  3.0d0, 4.58257569495584d0, 4.582575694955841d0, 3.0d0, 1.d00 /)
+
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_6, (28) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=6
@@ -243,10 +275,22 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ]
   cart_to_sphe_6 ( 2,13) = 1.2151388809514737933d0
   cart_to_sphe_6 ( 7,13) = -1.9764235376052370825d0
   cart_to_sphe_6 (16,13) = 1.2151388809514737933d0
+
+  cart_to_sphe_norm_6 = (/ 1.0d0, 3.3166247903554003d0, 3.3166247903554003d0, &
+  5.744562646538029d0, 9.949874371066201d0, 5.744562646538029d0, &
+  6.797058187186571d0, 15.198684153570666d0, 15.198684153570664d0, &
+  6.797058187186572d0, 5.744562646538029d0, 15.198684153570666d0, &
+  19.621416870348583d0, 15.198684153570666d0, 5.744562646538029d0, &
+  3.3166247903554003d0, 9.949874371066201d0, 15.198684153570664d0, &
+  15.198684153570666d0, 9.9498743710662d0, 3.3166247903554003d0, 1.0d0, &
+  3.3166247903554003d0, 5.744562646538029d0, 6.797058187186572d0, &
+  5.744562646538029d0, 3.3166247903554003d0, 1.d00 /)
+
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_7, (36) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=7
@@ -355,10 +399,25 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ]
   cart_to_sphe_7 ( 7,15) = -2.4456993503903949804d0
   cart_to_sphe_7 (16,15) = 1.96875d0
   cart_to_sphe_7 (29,15) = -0.64725984928774934788d0
+
+  cart_to_sphe_norm_7 = (/ 1.0d0, 3.6055512754639896d0, 3.605551275463989d0, &
+  6.904105059069327d0, 11.958260743101398d0, 6.904105059069326d0, &
+  9.26282894152753d0, 20.712315177207984d0, 20.71231517720798d0, &
+  9.26282894152753d0, 9.26282894152753d0, 24.507141816213494d0, &
+  31.63858403911275d0, 24.507141816213494d0, 9.262828941527529d0, &
+  6.904105059069327d0, 20.712315177207984d0, 31.63858403911275d0, &
+  31.63858403911275d0, 20.71231517720798d0, 6.904105059069327d0, &
+  3.6055512754639896d0, 11.958260743101398d0, 20.71231517720798d0, &
+  24.507141816213494d0, 20.71231517720798d0, 11.958260743101398d0, &
+  3.6055512754639896d0, 1.0d0, 3.605551275463989d0, 6.904105059069326d0, &
+  9.26282894152753d0, 9.262828941527529d0, 6.904105059069327d0, &
+  3.6055512754639896d0, 1.d00 /)
+
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_8, (45) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=8
@@ -506,10 +565,28 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ]
   cart_to_sphe_8 ( 7,17) = -2.9348392204684739765d0
   cart_to_sphe_8 (16,17) = 2.9348392204684739765d0
   cart_to_sphe_8 (29,17) = -1.2945196985754986958d0
+
+  cart_to_sphe_norm_8 = (/ 1.0d0, 3.872983346207417d0, 3.872983346207417d0, &
+  8.062257748298551d0, 13.964240043768942d0, 8.06225774829855d0, &
+  11.958260743101398d0, 26.739483914241877d0, 26.739483914241877d0, &
+  11.958260743101398d0, 13.55939315961975d0, 35.874782229304195d0, &
+  46.31414470763765d0, 35.874782229304195d0, 13.55939315961975d0, &
+  11.958260743101398d0, 35.874782229304195d0, 54.79963503528103d0, &
+  54.79963503528103d0, 35.874782229304195d0, 11.958260743101398d0, &
+  8.062257748298551d0, 26.739483914241877d0, 46.31414470763765d0, &
+  54.79963503528103d0, 46.314144707637645d0, 26.739483914241877d0, &
+  8.06225774829855d0, 3.872983346207417d0, 13.964240043768942d0, &
+  26.739483914241877d0, 35.874782229304195d0, 35.874782229304195d0, &
+  26.739483914241877d0, 13.96424004376894d0, 3.8729833462074166d0, 1.0d0, &
+  3.872983346207417d0, 8.06225774829855d0, 11.958260743101398d0, &
+  13.55939315961975d0, 11.958260743101398d0, 8.06225774829855d0, &
+  3.8729833462074166d0, 1.d0 /)
+
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ]
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_9, (55) ]
   implicit none
   BEGIN_DOC
 ! Spherical -> Cartesian Transformation matrix for l=9
@@ -703,5 +780,28 @@ BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ]
   cart_to_sphe_9 (16,19) = 4.1179360680974030877d0
   cart_to_sphe_9 (29,19) = -2.3781845426185916576d0
   cart_to_sphe_9 (46,19) = 0.60904939217552380708d0
+
+  cart_to_sphe_norm_9 = (/ 1.0d0, 4.1231056256176615d0, 4.1231056256176615d0, &
+  9.219544457292889d0, 15.968719422671313d0, 9.219544457292889d0, &
+  14.86606874731851d0, 33.24154027718933d0, 33.24154027718933d0, &
+  14.866068747318508d0, 18.635603405463275d0, 49.30517214248421d0, &
+  63.652703529910404d0, 49.30517214248421d0, 18.635603405463275d0, &
+  18.635603405463275d0, 55.90681021638982d0, 85.39906322671229d0, &
+  85.39906322671229d0, 55.90681021638983d0, 18.635603405463275d0, &
+  14.86606874731851d0, 49.30517214248421d0, 85.39906322671229d0, &
+  101.04553429023969d0, 85.3990632267123d0, 49.30517214248421d0, &
+  14.866068747318508d0, 9.219544457292889d0, 33.24154027718933d0, &
+  63.652703529910404d0, 85.39906322671229d0, 85.3990632267123d0, &
+  63.65270352991039d0, 33.24154027718933d0, 9.219544457292887d0, &
+  4.1231056256176615d0, 15.968719422671313d0, 33.24154027718933d0, &
+  49.30517214248421d0, 55.90681021638983d0, 49.30517214248421d0, &
+  33.24154027718933d0, 15.968719422671313d0, 4.1231056256176615d0, 1.0d0, &
+  4.1231056256176615d0, 9.219544457292889d0, 14.866068747318508d0, &
+  18.635603405463275d0, 18.635603405463275d0, 14.866068747318508d0, &
+  9.219544457292887d0, 4.1231056256176615d0, 1.d0 /)
+
 END_PROVIDER
 
+
+
+
diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f
index 2c6b8e7e..8b0edfbc 100644
--- a/src/ao_one_e_ints/ao_one_e_ints.irp.f
+++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f
@@ -45,3 +45,19 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)]
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals,(ao_sphe_num,ao_sphe_num)]
+&BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_diag,(ao_sphe_num)]
+  implicit none
+  integer :: i,j,n,l
+  BEGIN_DOC
+ ! One-electron Hamiltonian in the spherical |AO| basis.
+  END_DOC
+
+  ao_sphe_one_e_integrals = ao_sphe_integrals_n_e + ao_sphe_kinetic_integrals
+
+  do j = 1, ao_num
+    ao_sphe_one_e_integrals_diag(j) = ao_sphe_one_e_integrals(j,j)
+  enddo
+
+END_PROVIDER
+
diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f
index 747d9183..4e4d33f0 100644
--- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f
+++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f
@@ -1,9 +1,12 @@
  BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)]
+&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)]
 &BEGIN_PROVIDER [ integer, ao_sphe_num ]
   implicit none
   BEGIN_DOC
 ! Coefficients to go from cartesian to spherical coordinates in the current
 ! basis set
+!
+!  S_cart^-1 <cart|sphe>
   END_DOC
   integer :: i
   integer, external              :: ao_power_index
@@ -11,6 +14,7 @@
   integer                        :: prev
   prev = 0
   ao_cart_to_sphe_coef(:,:) = 0.d0
+  ao_cart_to_sphe_normalization(:) = 1.d0
   ! Assume order provided by ao_power_index
   i = 1
   ao_sphe_num = 0
@@ -19,6 +23,7 @@
       case (0)
         ao_sphe_num += 1
         ao_cart_to_sphe_coef(i,ao_sphe_num) = 1.d0
+        ao_cart_to_sphe_normalization(i) = 1.d0
         i += 1
       BEGIN_TEMPLATE
       case ($SHELL)
@@ -28,6 +33,9 @@
               ao_cart_to_sphe_coef(i+j-1,ao_sphe_num+k) = cart_to_sphe_$SHELL(j,k)
             enddo
           enddo
+          do j=1,size(cart_to_sphe_$SHELL,1)
+            ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j)
+          enddo
           i += size(cart_to_sphe_$SHELL,1)
           ao_sphe_num += size(cart_to_sphe_$SHELL,2)
         endif
@@ -47,28 +55,7 @@
     end select
   enddo
 
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, ao_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
- implicit none
- BEGIN_DOC
- ! |AO| overlap matrix in the spherical basis set
- END_DOC
- double precision, allocatable :: tmp(:,:)
- allocate (tmp(ao_sphe_num,ao_num))
-
- call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
-   ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
-   S_inv,size(ao_overlap,1), 0.d0, &
-   tmp, size(tmp,1))
-
- call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
-   tmp, size(tmp,1), &
-   ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, &
-   ao_sphe_overlap,size(ao_sphe_overlap,1))
-
- deallocate(tmp)
-
+print *, ao_cart_to_sphe_normalization(:)
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ]
@@ -77,9 +64,26 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ]
  ! Inverse of :c:data:`ao_cart_to_sphe_coef`
  END_DOC
 
- call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),&
-   ao_num,ao_sphe_num, &
-   ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1), lin_dep_cutoff)
+  ! Normalize
+  integer :: m,k
+  double precision, allocatable :: S(:,:), R(:,:), Rinv(:,:), Sinv(:,:)
+
+  k = size(ao_cart_to_sphe_coef,1)
+  m = size(ao_cart_to_sphe_coef,2)
+
+  allocate(S(k,k), R(k,m), Rinv(m,k), Sinv(k,k))
+
+  R(:,:) = ao_cart_to_sphe_coef(:,:)
+
+  call dgemm('N','T', m, m, k, 1.d0, R, k, R, k, 0.d0, S, m)
+  call get_pseudo_inverse(S, k, k, m, Sinv, k, 1.d-20)
+  call dgemm('T','T', m, m, k, 1.d0, R, k, Sinv, k, 0.d0, Rinv, m)
+
+  integer :: i
+  do i=1,ao_num
+    ao_cart_to_sphe_inv(:,i) = Rinv(:,i)  !/ ao_cart_to_sphe_normalization(i)
+  enddo
+
 END_PROVIDER
 
 
diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f
index e5341f6a..cc676c4b 100644
--- a/src/ao_one_e_ints/ao_overlap.irp.f
+++ b/src/ao_one_e_ints/ao_overlap.irp.f
@@ -308,3 +308,26 @@ BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num)  ]
 
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, ao_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! |AO| overlap matrix in the spherical basis set
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
+   ao_overlap,size(ao_overlap,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   tmp, size(tmp,1), &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
+   ao_sphe_overlap,size(ao_sphe_overlap,1))
+
+ deallocate(tmp)
+
+END_PROVIDER
+
diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f
index 49eb56ad..f450721d 100644
--- a/src/ao_one_e_ints/kin_ao_ints.irp.f
+++ b/src/ao_one_e_ints/kin_ao_ints.irp.f
@@ -190,3 +190,25 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)]
   endif
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, ao_sphe_kinetic_integrals, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! |AO| kinetic inntegrals matrix in the spherical basis set
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
+   ao_kinetic_integrals,size(ao_kinetic_integrals,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   tmp, size(tmp,1), &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
+   ao_sphe_kinetic_integrals,size(ao_sphe_kinetic_integrals,1))
+
+ deallocate(tmp)
+
+END_PROVIDER
diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f
index 776b5ec0..c1544a5d 100644
--- a/src/ao_one_e_ints/pot_ao_ints.irp.f
+++ b/src/ao_one_e_ints/pot_ao_ints.irp.f
@@ -609,3 +609,25 @@ double precision function V_r(n,alpha)
 end
 
 
+
+BEGIN_PROVIDER [ double precision, ao_sphe_integrals_n_e, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! |AO| VneVne  inntegrals matrix in the spherical basis set
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
+   ao_integrals_n_e,size(ao_integrals_n_e,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   tmp, size(tmp,1), &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
+   ao_sphe_integrals_n_e,size(ao_sphe_integrals_n_e,1))
+
+ deallocate(tmp)
+
+END_PROVIDER
diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
index e75ca056..851f26d8 100644
--- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
+++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
@@ -296,3 +296,67 @@ END_PROVIDER
  enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_local, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! |AO| pseudo_integrals_local matrix in the spherical basis set
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
+   ao_pseudo_integrals_local,size(ao_pseudo_integrals_local,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   tmp, size(tmp,1), &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
+   ao_sphe_pseudo_integrals_local,size(ao_sphe_pseudo_integrals_local,1))
+
+ deallocate(tmp)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_non_local, (ao_sphe_num,ao_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! |AO| pseudo_integrals_non_local matrix in the spherical basis set
+ END_DOC
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
+   ao_pseudo_integrals_non_local,size(ao_pseudo_integrals_non_local,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
+   tmp, size(tmp,1), &
+   ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
+   ao_sphe_pseudo_integrals_non_local,size(ao_sphe_pseudo_integrals_non_local,1))
+
+ deallocate(tmp)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals, (ao_sphe_num,ao_sphe_num)]
+  implicit none
+  BEGIN_DOC
+  ! Pseudo-potential integrals in the |AO| basis set.
+  END_DOC
+
+  ao_sphe_pseudo_integrals = 0.d0
+  if (do_pseudo) then
+    if (pseudo_klocmax > 0) then
+      ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_local
+    endif
+    if (pseudo_kmax > 0) then
+      ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_non_local
+    endif
+  endif
+
+END_PROVIDER
+
diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f
index a8b5441d..1eecca6c 100644
--- a/src/mo_basis/mos.irp.f
+++ b/src/mo_basis/mos.irp.f
@@ -347,11 +347,19 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
 end
 
 
-BEGIN_PROVIDER [ double precision, mo_coef_spherical]
+BEGIN_PROVIDER [ double precision, mo_sphe_coef, (ao_sphe_num, mo_num) ]
  implicit none
  BEGIN_DOC
  ! MO coefficients in the basis of spherical harmonics AOs.
  END_DOC
- 
+ double precision, allocatable :: tmp(:,:)
+ allocate (tmp(ao_sphe_num,ao_num))
+
+ call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_coef,ao_num, &
+   mo_coef,size(mo_coef,1), 0.d0, &
+   mo_sphe_coef, size(mo_sphe_coef,1))
+
+ deallocate (tmp)
 END_PROVIDER
 
diff --git a/src/trexio/EZFIO.cfg b/src/trexio/EZFIO.cfg
index 88828520..22dbea00 100644
--- a/src/trexio/EZFIO.cfg
+++ b/src/trexio/EZFIO.cfg
@@ -22,6 +22,13 @@ doc: If True, export MO coefficients
 interface: ezfio, ocaml, provider
 default: True
 
+[export_cartesian]
+type: logical
+doc: If False, export everything in the spherical AO basis
+interface: ezfio, ocaml, provider
+default: True
+
+
 [export_ao_one_e_ints]
 type: logical
 doc: If True, export one-electron integrals in AO basis
diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f
index cf1327b6..53b21dc9 100644
--- a/src/trexio/export_trexio_routines.irp.f
+++ b/src/trexio/export_trexio_routines.irp.f
@@ -77,6 +77,7 @@ subroutine export_trexio(update,full_path)
     rc = trexio_read_metadata_code_num(f(k), code_num)
     if (rc == TREXIO_ATTR_MISSING) then
       i = 1
+      code_num = 0
       code(:) = ""
     else
       rc = trexio_read_metadata_code(f(k), code, 64)
@@ -97,6 +98,7 @@ subroutine export_trexio(update,full_path)
     rc = trexio_read_metadata_author_num(f(k), author_num)
     if (rc == TREXIO_ATTR_MISSING) then
       i = 1
+      author_num = 0
       author(:) = ""
     else
       rc = trexio_read_metadata_author(f(k), author, 64)
@@ -305,35 +307,46 @@ subroutine export_trexio(update,full_path)
 
     print *, 'AOs'
 
-    rc = trexio_write_ao_num(f(1), ao_num)
-    call trexio_assert(rc, TREXIO_SUCCESS)
+    if (export_cartesian) then
+      rc = trexio_write_ao_cartesian(f(1), 1)
+      call trexio_assert(rc, TREXIO_SUCCESS)
 
-    rc = trexio_write_ao_cartesian(f(1), 1)
-    call trexio_assert(rc, TREXIO_SUCCESS)
+      rc = trexio_write_ao_num(f(1), ao_num)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_shell(f(1), ao_shell)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      if (ezfio_convention >= 20250211) then
+        rc = trexio_write_ao_normalization(f(1), ao_coef_normalization_factor)
+      else
+
+        allocate(factor(ao_num))
+        do i=1,ao_num
+          l = ao_first_of_shell(ao_shell(i))
+          factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
+        enddo
+        rc = trexio_write_ao_normalization(f(1), factor)
+        deallocate(factor)
+      endif
+
+      call trexio_assert(rc, TREXIO_SUCCESS)
 
-    rc = trexio_write_ao_shell(f(1), ao_shell)
-    call trexio_assert(rc, TREXIO_SUCCESS)
 
-    if (ezfio_convention >= 20250211) then
-      rc = trexio_write_ao_normalization(f(1), ao_coef_normalization_factor)
     else
-      integer :: pow0(3), powA(3), nz
-      double precision :: normA, norm0, C_A(3), overlap_x, overlap_z, overlap_y, c
-      nz=100
+      rc = trexio_write_ao_cartesian(f(1), 0)
+      call trexio_assert(rc, TREXIO_SUCCESS)
 
-      C_A(1) = 0.d0
-      C_A(2) = 0.d0
-      C_A(3) = 0.d0
+      rc = trexio_write_ao_num(f(1), ao_sphe_num)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_shell(f(1), ao_sphe_shell)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_normalization(f(1), ao_sphe_coef_normalization_factor)
+      call trexio_assert(rc, TREXIO_SUCCESS)
 
-      allocate(factor(ao_num))
-      do i=1,ao_num
-        l = ao_first_of_shell(ao_shell(i))
-        factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
-      enddo
-      rc = trexio_write_ao_normalization(f(1), factor)
-      deallocate(factor)
     endif
-    call trexio_assert(rc, TREXIO_SUCCESS)
 
   endif
 
@@ -341,23 +354,45 @@ subroutine export_trexio(update,full_path)
 ! ------------------
 
   if (export_ao_one_e_ints) then
-    print *, 'AO one-e integrals'
 
-    rc = trexio_write_ao_1e_int_overlap(f(1),ao_overlap)
-    call trexio_assert(rc, TREXIO_SUCCESS)
+    double precision, allocatable :: tmp_ao(:,:)
+    if (export_cartesian) then
+      print *, 'AO one-e integrals (cartesian)'
 
-    rc = trexio_write_ao_1e_int_kinetic(f(1),ao_kinetic_integrals)
-    call trexio_assert(rc, TREXIO_SUCCESS)
-
-    rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_integrals_n_e)
-    call trexio_assert(rc, TREXIO_SUCCESS)
-
-    if (do_pseudo) then
-      rc = trexio_write_ao_1e_int_ecp(f(1), ao_pseudo_integrals_local + ao_pseudo_integrals_non_local)
+      rc = trexio_write_ao_1e_int_overlap(f(1),ao_overlap)
       call trexio_assert(rc, TREXIO_SUCCESS)
-    endif
 
-    rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_one_e_integrals)
+      rc = trexio_write_ao_1e_int_kinetic(f(1),ao_kinetic_integrals)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_integrals_n_e)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      if (do_pseudo) then
+        rc = trexio_write_ao_1e_int_ecp(f(1), ao_pseudo_integrals)
+        call trexio_assert(rc, TREXIO_SUCCESS)
+      endif
+
+      rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_one_e_integrals)
+    else
+      print *, 'AO one-e integrals (spherical)'
+
+      rc = trexio_write_ao_1e_int_overlap(f(1),ao_sphe_overlap)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_1e_int_kinetic(f(1),ao_sphe_kinetic_integrals)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      rc = trexio_write_ao_1e_int_potential_n_e(f(1),ao_sphe_integrals_n_e)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+
+      if (do_pseudo) then
+        rc = trexio_write_ao_1e_int_ecp(f(1), ao_sphe_pseudo_integrals)
+        call trexio_assert(rc, TREXIO_SUCCESS)
+      endif
+
+      rc = trexio_write_ao_1e_int_core_hamiltonian(f(1),ao_sphe_one_e_integrals)
+    endif
     call trexio_assert(rc, TREXIO_SUCCESS)
   end if
 
@@ -465,8 +500,13 @@ subroutine export_trexio(update,full_path)
       call trexio_assert(rc, TREXIO_SUCCESS)
     enddo
 
-    rc = trexio_write_mo_coefficient(f(1), mo_coef)
-    call trexio_assert(rc, TREXIO_SUCCESS)
+    if (export_cartesian) then
+      rc = trexio_write_mo_coefficient(f(1), mo_coef)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    else
+      rc = trexio_write_mo_coefficient(f(1), mo_sphe_coef)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
 
     if ( (trim(mo_label) == 'Canonical').and. &
          (export_mo_two_e_ints_cholesky.or.export_mo_two_e_ints) ) then
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index 4e7ca87d..629a998b 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1392,15 +1392,6 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff)
 
   call dgemm('T', 'T', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1))
 
-!  C = 0.d0
-!  do i=1,m
-!    do j=1,n
-!      do k=1,n_svd
-!        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
-!      enddo
-!    enddo
-!  enddo
-
   deallocate(U,D,Vt,work,A_tmp)
 
 end