From 7a0194d57690e67e44218aca4a75a1a3e3c4951c Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 27 Nov 2024 14:19:53 +0100
Subject: [PATCH] Fixed problems with trexio files coming from cp2k

---
 scripts/qp_import_trexio.py | 13 ++++++++-----
 src/ao_basis/aos.irp.f      | 24 ++++++++++++------------
 2 files changed, 20 insertions(+), 17 deletions(-)

diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py
index 70551c8f..c25a27a9 100755
--- a/scripts/qp_import_trexio.py
+++ b/scripts/qp_import_trexio.py
@@ -160,7 +160,6 @@ def write_ezfio(trexio_filename, filename):
             ezfio.set_basis_shell_ang_mom(ang_mom)
             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)
 
             nucl_shell_num = []
             prev = None
@@ -194,6 +193,10 @@ def write_ezfio(trexio_filename, filename):
 
             shell_factor = trexio.read_basis_shell_factor(trexio_file)
             prim_factor  = trexio.read_basis_prim_factor(trexio_file)
+            for i,p in enumerate(prim_factor):
+                coefficient[i] *= prim_factor[i]
+            ezfio.set_ao_basis_primitives_normalized(False)
+            ezfio.set_basis_prim_coef(coefficient)
 
         elif basis_type.lower() == "numerical":
 
@@ -245,7 +248,6 @@ def write_ezfio(trexio_filename, filename):
             ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
 
             shell_factor = trexio.read_basis_shell_factor(trexio_file)
-            prim_factor  = [1.]*prim_num
         else:
            raise TypeError
 
@@ -387,10 +389,11 @@ def write_ezfio(trexio_filename, filename):
 
       # Renormalize MO coefs if needed
       if trexio.has_ao_normalization(trexio_file_cart):
+        ezfio.set_ao_basis_ao_normalized(False)
         norm = trexio.read_ao_normalization(trexio_file_cart)
-        for j in range(mo_num):
-          for i,f in enumerate(norm):
-             MoMatrix[i,j] *= f
+#        for j in range(mo_num):
+#           for i,f in enumerate(norm):
+#              MoMatrix[i,j] *= f
       ezfio.set_mo_basis_mo_coef(MoMatrix)
 
       mo_occ = [ 0. for i in range(mo_num) ]
diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f
index c9213f47..34853398 100644
--- a/src/ao_basis/aos.irp.f
+++ b/src/ao_basis/aos.irp.f
@@ -55,9 +55,6 @@ END_PROVIDER
 
   do i=1,ao_num
 
-!    powA(1) = ao_power(i,1) +  ao_power(i,2) +  ao_power(i,3)
-!    powA(2) = 0
-!    powA(3) = 0
     powA(1) = ao_power(i,1)
     powA(2) = ao_power(i,2)
     powA(3) = ao_power(i,3)
@@ -76,16 +73,19 @@ END_PROVIDER
     endif
 
     ! Normalization of the contracted basis functions
-    if (ao_normalized) then
-      norm = 0.d0
-      do j=1,ao_prim_num(i)
-        do k=1,ao_prim_num(i)
-          call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
-          norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
-        enddo
+    norm = 0.d0
+    do j=1,ao_prim_num(i)
+      do k=1,ao_prim_num(i)
+        call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+        norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
+      enddo
+    enddo
+    ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
+
+    if (.not.ao_normalized) then
+      do j=1,ao_prim_num(i)
+        ao_coef_normalized(i,j) = ao_coef_normalized(i,j) * ao_coef_normalization_factor(i)
       enddo
-      ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
-    else
       ao_coef_normalization_factor(i) = 1.d0
     endif
   enddo