Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

2025-04-25 17:54:44 +02:00 · 2025-03-20 15:35:40 +01:00 · 2025-03-20 15:35:40 +01:00 · 7a2e24765c
commit 7a2e24765c
parent b1fa902513 66ce7bba22
19 changed files with 582 additions and 132 deletions
--- a/etc/openmp.rc
+++ b/etc/openmp.rc
@ -1 +1 @@
-export OMP_NESTED=True
+#export OMP_NESTED=True
--- a/external/irpf90
+++ b/external/irpf90
@ -1 +1 @@
-Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
+Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b
--- a/plugins/local/basis_correction/weak_corr_func.irp.f
+++ b/plugins/local/basis_correction/weak_corr_func.irp.f
@ -1,45 +1,73 @@

- BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)]                                                                      
+ BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)]
 BEGIN_DOC
-! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) , 
+! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) ,
 !
 ! see equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714
 END_DOC
 implicit none
 integer :: ipoint,istate
- double precision :: rho_a, rho_b, ec
+ double precision :: rho_a, rho_b, ec, rho, p2
 double precision :: wall0,wall1,weight,mu
 logical :: dospin
+ double precision, external :: g0_UEG_mu_inf
 dospin = .true. ! JT dospin have to be set to true for open shell
 print*,'Providing ecmd_lda_mu_of_r ...'

 ecmd_lda_mu_of_r = 0.d0
 call wall_time(wall0)
- do istate = 1, N_states
-  do ipoint = 1, n_points_final_grid
-   ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
-   mu     = mu_of_r_prov(ipoint,istate)
-   weight = final_weight_at_r_vector(ipoint)
-   rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
-   rho_b =  one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
-   ! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
-   call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
-   if(isnan(ec))then
-    print*,'ec is nan'
-    stop
-   endif
-   ecmd_lda_mu_of_r(istate) += weight * ec
+ if (mu_of_r_potential.EQ."proj_DUMMY")then
+  do istate = 1, N_states
+    do ipoint = 1, n_points_final_grid
+    ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
+    mu     = mu_of_r_prov(ipoint,istate)
+    weight = final_weight_at_r_vector(ipoint)
+    rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
+    rho_b =  one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
+    rho = rho_a + rho_b
+!    p2 = rho_a*rho_b
+!  We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
+    p2 = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
+    rho_a = 0.5d0*(rho + dsqrt(rho*rho - 4.d0*p2))
+    rho_b = 0.5d0*(rho - dsqrt(rho*rho - 4.d0*p2))
+!    rho_a = 0.5d0*rho
+!    rho_b = 0.5d0*rho
+    ! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
+    call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
+    if(isnan(ec))then
+      print*,'ec is nan'
+      stop
+    endif
+    ecmd_lda_mu_of_r(istate) += weight * ec
+    enddo
  enddo
- enddo
+ else
+  do istate = 1, N_states
+    do ipoint = 1, n_points_final_grid
+    ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
+    mu     = mu_of_r_prov(ipoint,istate)
+    weight = final_weight_at_r_vector(ipoint)
+    rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
+    rho_b =  one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
+    ! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
+    call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
+    if(isnan(ec))then
+      print*,'ec is nan'
+      stop
+    endif
+    ecmd_lda_mu_of_r(istate) += weight * ec
+    enddo
+  enddo
+ endif
 call wall_time(wall1)
 print*,'Time for ecmd_lda_mu_of_r :',wall1-wall0
- END_PROVIDER 
+END_PROVIDER


 BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
 BEGIN_DOC
 ! ecmd_pbe_ueg_mu_of_r = multi-determinantal Ecmd within the PBE-UEG approximation with mu(r) ,
- ! 
+ !
 ! see Eqs. 13-14b in Phys.Chem.Lett.2019, 10, 2931   2937;  https://pubs.acs.org/doi/10.1021/acs.jpclett.9b01176
 !
 ! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638)
@ -50,10 +78,10 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
 double precision :: weight
 integer :: ipoint,istate
 double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
- double precision :: g0_UEG_mu_inf
+ double precision, external :: g0_UEG_mu_inf

 ecmd_pbe_ueg_mu_of_r = 0.d0
-  
+
 print*,'Providing ecmd_pbe_ueg_mu_of_r ...'
 call wall_time(wall0)
 do istate = 1, N_states
@ -69,7 +97,7 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
   grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)

 !  We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
-   on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)   
+   on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)

 ! The form of interpolated (mu=0 ---> mu=infinity) functional originally introduced in JCP, 150, 084103 1-10 (2019)
   call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE)
--- a/scripts/qp_import_trexio.py
+++ b/scripts/qp_import_trexio.py
@ -210,6 +210,7 @@ def write_ezfio(trexio_filename, filename):
            nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
            exponent    = [1.]*prim_num
            coefficient = [1.]*prim_num
+            prim_factor = [1.]*prim_num
            shell_index = [i for i in range(shell_num)]
            ao_shell    = trexio.read_ao_shell(trexio_file)

@ -221,6 +222,9 @@ def write_ezfio(trexio_filename, filename):
            ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
            ezfio.set_basis_prim_expo(exponent)
            ezfio.set_basis_prim_coef(coefficient)
+            ezfio.set_basis_prim_normalization_factor(prim_factor)
+            ezfio.set_basis_primitives_normalized(True)
+            ezfio.set_basis_ao_normalized(False)

            nucl_shell_num = []
            prev = None
@ -271,12 +275,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

@ -284,7 +287,7 @@ def write_ezfio(trexio_filename, filename):
    ezfio.set_ao_basis_ao_num(ao_num)


-    if cartesian and basis_type.lower() == "gaussian" and shell_num > 0:
+    if cartesian and basis_type.lower() in ["gaussian", "numerical"] and shell_num > 0:
        ao_shell    = trexio.read_ao_shell(trexio_file_cart)
        at = [ nucl_index[i]+1 for i in ao_shell ]
        ezfio.set_ao_basis_ao_nucl(at)
@ -521,11 +524,11 @@ def write_ezfio(trexio_filename, filename):
    alpha = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(alpha_s) ][::-1]
    beta  = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(beta_s ) ][::-1]
    ezfio.set_determinants_bit_kind(8)
-    ezfio.set_determinants_n_int(1+mo_num//64)
+    ezfio.set_determinants_n_int(1+(mo_num-1)//64)
    ezfio.set_determinants_n_det(1)
    ezfio.set_determinants_n_states(1)
-    ezfio.set_determinants_psi_det(alpha+beta)
    ezfio.set_determinants_psi_coef([[1.0]])
+    ezfio.set_determinants_psi_det(alpha+beta)
    print("OK")


--- a/src/ao_basis/aos.irp.f
+++ b/src/ao_basis/aos.irp.f
@ -20,7 +20,35 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ]
     ao_shell(k) = i
   enddo
 enddo
+END_PROVIDER

+BEGIN_PROVIDER [ integer, ao_sphe_num ]
+ implicit none
+ BEGIN_DOC
+ ! Number of spherical AOs
+ END_DOC
+ integer :: n, i
+ ao_sphe_num=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)
+   ao_sphe_num += 2*n+1
+ 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 +161,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) ]
--- 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

+
+
+
--- 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
+
--- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f
+++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f
@ -1,35 +1,42 @@
 BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)]
-&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ]
+&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_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
  integer                        :: ibegin,j,k
-  integer                        :: prev
+  integer                        :: prev, ao_sphe_count
  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_cart_to_sphe_num = 0
+  ao_sphe_count = 0
  do while (i <= ao_num)
    select case ( ao_l(i) )
      case (0)
-        ao_cart_to_sphe_num += 1
-        ao_cart_to_sphe_coef(i,ao_cart_to_sphe_num) = 1.d0
+        ao_sphe_count += 1
+        ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0
+        ao_cart_to_sphe_normalization(i) = 1.d0
        i += 1
      BEGIN_TEMPLATE
      case ($SHELL)
        if (ao_power(i,1) == $SHELL) then
          do k=1,size(cart_to_sphe_$SHELL,2)
            do j=1,size(cart_to_sphe_$SHELL,1)
-              ao_cart_to_sphe_coef(i+j-1,ao_cart_to_sphe_num+k) = cart_to_sphe_$SHELL(j,k)
+              ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+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_cart_to_sphe_num += size(cart_to_sphe_$SHELL,2)
+          ao_sphe_count += size(cart_to_sphe_$SHELL,2)
        endif
      SUBST [ SHELL ]
        1;;
@ -47,22 +54,25 @@
    end select
  enddo

+  if (ao_sphe_count /= ao_sphe_num) then
+    call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num")
+  endif
 END_PROVIDER

-BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ]
+BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
 implicit none
 BEGIN_DOC
- ! |AO| overlap matrix in the spherical basis set
+ ! T^T . S . T
 END_DOC
 double precision, allocatable :: S(:,:)
- allocate (S(ao_cart_to_sphe_num,ao_num))
+ allocate (S(ao_sphe_num,ao_num))

- call dgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, 1.d0, &
+ 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), &
   ao_overlap,size(ao_overlap,1), 0.d0, &
   S, size(S,1))

- call dgemm('N','N',ao_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, 1.d0, &
+ call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
   S, size(S,1), &
   ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, &
   ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1))
@ -71,15 +81,33 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num

 END_PROVIDER

-BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_num) ]
+BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ]
 implicit none
 BEGIN_DOC
 ! 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_cart_to_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-12)
+  call dgemm('T','N', 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)
+  enddo
+
 END_PROVIDER


@ -120,17 +148,17 @@ END_PROVIDER

    double precision, allocatable :: S(:,:)

-    allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num))
+    allocate(S(ao_sphe_num,ao_sphe_num))
    S = 0.d0
-    do i=1,ao_cart_to_sphe_num
+    do i=1,ao_sphe_num
      S(i,i) = 1.d0
    enddo

-    ao_ortho_canonical_num = ao_cart_to_sphe_num
+    ao_ortho_canonical_num = ao_sphe_num
    call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
-      ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff)
+      ao_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff)

-    call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, &
+    call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_sphe_num, 1.d0, &
      ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), &
      S, size(S,1), &
      0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1))
@ -167,3 +195,4 @@ BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonica
    enddo
  enddo
 END_PROVIDER
+
--- 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('N','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','T',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
+
--- 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('N','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','T',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
--- 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('N','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','T',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
--- 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
+
--- a/src/ao_two_e_ints/two_e_integrals.irp.f
+++ b/src/ao_two_e_ints/two_e_integrals.irp.f
@ -649,11 +649,20 @@ double precision function general_primitive_integral(dim,            &
 !  call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
  if (ior(n_pt_tmp,n_Iz) >= 0) then
    ! Bottleneck here
-    do ib=0,n_pt_tmp
-      do ic = 0,n_Iz
-        d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
+    if (ic > ib) then
+      do ib=0,n_pt_tmp
+          d1(ib:) = d1(ib:) + Iz_pol(:) * d_poly(ib)
      enddo
-    enddo
+    else
+      do ic = 0,n_Iz
+        d1(ic:) = d1(ic:) + Iz_pol(ic) * d_poly(:)
+      enddo
+    endif
+!    do ib=0,n_pt_tmp
+!      do ic = 0,n_Iz
+!        d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
+!      enddo
+!    enddo

    do n_pt_out = n_pt_tmp+n_Iz, 0, -1
      if (d1(n_pt_out) /= 0.d0) exit
--- a/src/mo_basis/mos.irp.f
+++ b/src/mo_basis/mos.irp.f
@ -346,3 +346,37 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
  deallocate(T)
 end

+
+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, dimension (:,:) :: C0, S, tmp
+ allocate(C0(ao_sphe_num,mo_num), S(mo_num,mo_num), tmp(mo_num,ao_sphe_num))
+
+ call dgemm('T','N',ao_sphe_num,mo_num,ao_num, 1.d0, &
+   ao_cart_to_sphe_coef,ao_num, &
+   mo_coef,size(mo_coef,1), 0.d0, &
+   C0, size(C0,1))
+
+ ! C0^T S S^0
+ call dgemm('T','N',mo_num,ao_sphe_num,ao_sphe_num, 1.d0, &
+   C0,size(C0,1), &
+   ao_sphe_overlap,size(ao_sphe_overlap,1), 0.d0, &
+   tmp, size(tmp,1))
+
+ call dgemm('N','N',mo_num,mo_num,ao_sphe_num, 1.d0, &
+   tmp, size(tmp,1), &
+   C0,size(C0,1), 0.d0, &
+   S, size(S,1))
+
+ integer :: m
+ m = ao_sphe_num
+ mo_sphe_coef = C0
+ call ortho_lowdin(S,size(S,1),mo_num,mo_sphe_coef,ao_sphe_num,m,1.d-10)
+
+
+ deallocate (tmp, S, C0)
+END_PROVIDER
+
--- a/src/scf_utils/roothaan_hall_scf.irp.f
+++ b/src/scf_utils/roothaan_hall_scf.irp.f
@ -222,7 +222,7 @@ END_DOC
  endif


-  ! Identify degenerate MOs and force them to be on the axes
+  ! Identify degenerate MOs and combine them to force them to be on the axes
  allocate(S(ao_num,ao_num))
  i=1
  do while (i<mo_num)
--- 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
--- a/src/trexio/export_trexio_routines.irp.f
+++ b/src/trexio/export_trexio_routines.irp.f
@ -72,11 +72,12 @@ subroutine export_trexio(update,full_path)
  character*(64) :: code(100), author(100), user
  character*(64), parameter :: qp2_code = "QuantumPackage"

-  call getenv("USER",user)
+  call getenv('USER',user)
  do k=1,N_states
    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)
@ -94,24 +95,27 @@ subroutine export_trexio(update,full_path)
      call trexio_assert(rc, TREXIO_SUCCESS)
    endif

-    rc = trexio_read_metadata_author_num(f(k), author_num)
-    if (rc == TREXIO_ATTR_MISSING) then
-      i = 1
-      author(:) = ""
-    else
-      rc = trexio_read_metadata_author(f(k), author, 64)
-      do i=1, author_num
-        if (trim(author(i)) == trim(user)) then
-          exit
-        endif
-      enddo
-    endif
-    if (i == author_num+1) then
-      author(i) = user
-      rc = trexio_write_metadata_author_num(f(k), i)
-      call trexio_assert(rc, TREXIO_SUCCESS)
-      rc = trexio_write_metadata_author(f(k), author, 64)
-      call trexio_assert(rc, TREXIO_SUCCESS)
+    if (trim(user) /= '') then
+      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)
+        do i=1, author_num
+          if (trim(author(i)) == trim(user)) then
+            exit
+          endif
+        enddo
+      endif
+      if (i == author_num+1) then
+        author(i) = user
+        rc = trexio_write_metadata_author_num(f(k), i)
+        call trexio_assert(rc, TREXIO_SUCCESS)
+        rc = trexio_write_metadata_author(f(k), author, 64)
+        call trexio_assert(rc, TREXIO_SUCCESS)
+      endif
    endif
  enddo

@ -305,35 +309,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 +356,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 +502,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
--- a/src/trexio/import_trexio_integrals.irp.f
+++ b/src/trexio/import_trexio_integrals.irp.f
@ -119,6 +119,28 @@ subroutine run(f)
    call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read')
  endif

+  ! Some codes only provide ao_1e_int_core_hamiltonian rather than
+  ! kinetic and nuclear potentials separately, so we need to work
+  ! around that. This is needed for non-GTO basis sets since some QP
+  ! functions will try to calculate these matrices from the nonexisting 
+  ! GTO basis if they are not set.
+  if (trexio_has_ao_1e_int_core_hamiltonian(f) == TREXIO_SUCCESS .and. &
+      trexio_has_ao_1e_int_potential_n_e(f) /= TREXIO_SUCCESS .and. &
+      trexio_has_ao_1e_int_kinetic(f) /= TREXIO_SUCCESS) then
+    rc = trexio_read_ao_1e_int_core_hamiltonian(f, A)
+    if (rc /= TREXIO_SUCCESS) then
+      print *, irp_here
+      print *, 'Error reading AO core Hamiltonian.'
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      stop -1
+    endif
+    call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A)
+    call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read')
+    A=0.d0
+    call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A)
+    call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
+  endif
+
  deallocate(A,B)

  ! AO 2e integrals
--- 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