Merge pull request #363 from QuantumPackage/fix-trexio

Fixed TREXIO norm problems

ocaml/Input_ao_basis.ml

@@ -13,8 +13,6 @@ module Ao_basis : sig
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-      ao_normalized   : bool;
-      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
   val read : unit -> t option
@@ -36,8 +34,6 @@ end = struct
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
-      ao_normalized   : bool;
-      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
 
@@ -133,24 +129,6 @@ end = struct
     Ezfio.get_ao_basis_ao_cartesian ()
   ;;
 
-  let read_ao_normalized () =
-    if not (Ezfio.has_ao_basis_ao_normalized()) then
-       get_default "ao_normalized"
-       |> bool_of_string
-       |> Ezfio.set_ao_basis_ao_normalized
-    ;
-    Ezfio.get_ao_basis_ao_normalized ()
-  ;;
-
-  let read_primitives_normalized () =
-    if not (Ezfio.has_ao_basis_primitives_normalized()) then
-       get_default "primitives_normalized"
-       |> bool_of_string
-       |> Ezfio.set_ao_basis_primitives_normalized
-    ;
-    Ezfio.get_ao_basis_primitives_normalized ()
-  ;;
-
   let to_long_basis b =
     let ao_num = AO_number.to_int b.ao_num in
     let gto_array = Array.init (AO_number.to_int b.ao_num)
@@ -213,8 +191,6 @@ end = struct
          ao_coef         ;
          ao_expo         ;
          ao_cartesian    ;
-         ao_normalized   ;
-         primitives_normalized ;
        } = b
      in
      write_md5 b ;
@@ -247,8 +223,6 @@ end = struct
      ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
 
      Ezfio.set_ao_basis_ao_cartesian(ao_cartesian);
-     Ezfio.set_ao_basis_ao_normalized(ao_normalized);
-     Ezfio.set_ao_basis_primitives_normalized(primitives_normalized);
 
      let ao_coef =
       Array.to_list ao_coef
@@ -280,8 +254,6 @@ end = struct
             ao_coef         = read_ao_coef () ;
             ao_expo         = read_ao_expo () ;
             ao_cartesian    = read_ao_cartesian () ;
-            ao_normalized   = read_ao_normalized () ;
-            primitives_normalized   = read_primitives_normalized () ;
           }
         in
         to_md5 result
@@ -392,8 +364,6 @@ end = struct
       { ao_basis = name ;
         ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ;
         ao_power ; ao_coef ; ao_expo ; ao_cartesian ;
-        ao_normalized = bool_of_string @@ get_default "ao_normalized";
-        primitives_normalized = bool_of_string @@ get_default "primitives_normalized";
         }
   ;;
 
@@ -448,14 +418,6 @@ Cartesian coordinates (6d,10f,...) ::
 
   ao_cartesian = %s
 
-Use normalized primitive functions ::
-
-  primitives_normalized = %s
-
-Use normalized basis functions ::
-
-  ao_normalized = %s
-
 Basis set (read-only) ::
 
 %s
@@ -469,8 +431,6 @@ Basis set (read-only) ::
 
 "   (AO_basis_name.to_string b.ao_basis)
     (string_of_bool b.ao_cartesian)
-    (string_of_bool b.primitives_normalized)
-    (string_of_bool b.ao_normalized)
     (Basis.to_string short_basis
        |> String_ext.split ~on:'\n'
        |> list_map (fun x-> "  "^x)
@@ -507,8 +467,6 @@ ao_power                = %s
 ao_coef                 = %s
 ao_expo                 = %s
 ao_cartesian            = %s
-ao_normalized           = %s
-primitives_normalized   = %s
 md5                     = %s
 "
     (AO_basis_name.to_string b.ao_basis)
@@ -525,8 +483,6 @@ md5                     = %s
     (b.ao_expo  |> Array.to_list |> list_map AO_expo.to_string
       |> String.concat ", ")
     (b.ao_cartesian |> string_of_bool)
-    (b.ao_normalized |> string_of_bool)
-    (b.primitives_normalized |> string_of_bool)
     (to_md5 b |> MD5.to_string )
 
   ;;

scripts/qp_import_trexio.py

@@ -191,11 +191,14 @@ 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  = trexio.read_basis_prim_factor(trexio_file)
-            for i,p in enumerate(prim_factor):
+            ezfio.set_basis_prim_normalization_factor(prim_factor)
-                coefficient[i] *= prim_factor[i]
+            ezfio.set_basis_primitives_normalized(True)
-            ezfio.set_ao_basis_primitives_normalized(False)
+            ezfio.set_basis_ao_normalized(False)
+
+            shell_factor = trexio.read_basis_shell_factor(trexio_file)
+            for i, shell_idx in enumerate(shell_index):
+               coefficient[i] *= shell_factor[shell_idx]
             ezfio.set_basis_prim_coef(coefficient)
 
         elif basis_type.lower() == "numerical":
@@ -247,7 +250,6 @@ def write_ezfio(trexio_filename, filename):
             ezfio.set_basis_shell_index([x+1 for x in shell_index])
             ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
 
-            shell_factor = trexio.read_basis_shell_factor(trexio_file)
         else:
            raise TypeError
 
@@ -317,13 +319,18 @@ def write_ezfio(trexio_filename, filename):
                 exponent.append(expo[i])
                 num_prim.append(num_prim0[i])
 
-        print (len(coefficient), ao_num)
         assert (len(coefficient) == ao_num)
+
         ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
         ezfio.set_ao_basis_ao_prim_num(num_prim)
 
         prim_num_max = max( [ len(x) for x in coefficient ] )
 
+        ao_normalization = trexio.read_ao_normalization(trexio_file_cart)
+        for i, coef in enumerate(coefficient):
+           for j in range(len(coef)):
+              coef[j] *= ao_normalization[i]
+
         for i in range(ao_num):
             coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)]
             exponent   [i] += [0. for j in range(len(exponent[i]), prim_num_max)]
@@ -338,7 +345,6 @@ def write_ezfio(trexio_filename, filename):
                 coef.append(coefficient[j])
                 expo.append(exponent[j])
 
-#        ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
         ezfio.set_ao_basis_ao_coef(coef)
         ezfio.set_ao_basis_ao_expo(expo)
 
@@ -388,14 +394,6 @@ def write_ezfio(trexio_filename, filename):
 
       # Read coefs from temporary cartesian file created in the AO section
       MoMatrix = trexio.read_mo_coefficient(trexio_file_cart)
-
-      # 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
       ezfio.set_mo_basis_mo_coef(MoMatrix)
 
       mo_occ = [ 0. for i in range(mo_num) ]

src/ao_basis/EZFIO.cfg

@@ -55,18 +55,6 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...)
 interface: ezfio, provider
 default: false
 
-[ao_normalized]
-type: logical
-doc: Use normalized basis functions
-interface: ezfio, provider
-default: true
-
-[primitives_normalized]
-type: logical
-doc: Use normalized primitive functions
-interface: ezfio, provider
-default: true
-
 [use_cgtos]
 type: logical
 doc: If true, use cgtos for AO integrals

src/ao_basis/aos.irp.f

@@ -82,10 +82,11 @@ END_PROVIDER
     enddo
     ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
 
-    if (.not.ao_normalized) then
+    if (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
+    else
       ao_coef_normalization_factor(i) = 1.d0
     endif
   enddo

src/basis/EZFIO.cfg

@@ -73,3 +73,15 @@ size: (basis.prim_num)
 interface: ezfio, provider
 
 
+[primitives_normalized]
+type: logical
+doc: If true, assume primitive basis functions are normalized
+interface: ezfio, provider, ocaml
+default: true
+
+[ao_normalized]
+type: logical
+doc: If true, normalize the basis functions
+interface: ezfio, provider, ocaml
+default: false
+

src/basis/basis.irp.f

@@ -6,6 +6,11 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
 
   logical                        :: has
   PROVIDE ezfio_filename
+  if (.not.ao_normalized) then
+    shell_normalization_factor = 1.d0
+    return
+  endif
+
   if (mpi_master) then
     if (size(shell_normalization_factor) == 0) return
 
@@ -70,6 +75,12 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
 
   logical                        :: has
   PROVIDE ezfio_filename
+
+  if (.not.primitives_normalized) then
+    prim_normalization_factor(:) = 1.d0
+    return
+  endif
+
   if (mpi_master) then
     if (size(prim_normalization_factor) == 0) return
 
@@ -95,9 +106,9 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
 
         do k=1, prim_num
           if (shell_index(k) /= i) cycle
-            call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
+          call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
-              powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+            powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
-            prim_normalization_factor(k) = 1.d0/dsqrt(norm)
+          prim_normalization_factor(k) = 1.d0/dsqrt(norm)
         enddo
       enddo
 

src/mo_two_e_ints/map_integrals.irp.f

@@ -143,7 +143,7 @@ double precision function get_two_e_integral(i,j,k,l,map)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx
-  integer                        :: ii
+  integer                        :: ii, kk
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
 
@@ -179,8 +179,11 @@ double precision function get_two_e_integral(i,j,k,l,map)
 
       double precision, external :: ddot
       get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1)
-!       double precision, external :: get_from_mo_cholesky_cache
+
-!       get_two_e_integral = get_from_mo_cholesky_cache(i,j,k,l,.False.)
+!      get_two_e_integral = 0.d0
+!      do kk=1,cholesky_mo_num
+!        get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!      enddo
 
     else
 
@@ -516,11 +519,46 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
 
     if (do_mo_cholesky) then
 
-      double precision, external :: ddot
+      if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. &
-      do i=1,sze
+           (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then
-        out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+
-                                           cholesky_mo_transp(1,i,l), 1)
+        double precision, external :: ddot
-      enddo
+        integer :: kk
+
+        do i=1,mo_integrals_cache_min-1
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+        do i=mo_integrals_cache_min,mo_integrals_cache_max
+          out_val(i) = get_two_e_integral_cache(i,i,k,l)
+        enddo
+
+        do i=mo_integrals_cache_max, sze
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+      else
+
+        do i=1,sze
+          out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, &
+                                             cholesky_mo_transp(1,i,l), 1)
+!          out_val(i) = 0.d0
+!          do kk=1,cholesky_mo_num
+!            out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l)
+!          enddo
+        enddo
+
+      endif
 
     else
 

src/trexio/export_trexio_routines.irp.f

@@ -271,11 +271,7 @@ subroutine export_trexio(update,full_path)
      call trexio_assert(rc, TREXIO_SUCCESS)
 
      allocate(factor(shell_num))
-!     if (ao_normalized) then
+     factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-       factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-!     else
-!       factor(1:shell_num) = 1.d0
-!     endif
      rc = trexio_write_basis_shell_factor(f(1), factor)
      call trexio_assert(rc, TREXIO_SUCCESS)
 
@@ -291,11 +287,12 @@ subroutine export_trexio(update,full_path)
     call trexio_assert(rc, TREXIO_SUCCESS)
 
     allocate(factor(prim_num))
-!    if (primitives_normalized) then
+    if (primitives_normalized) then
       factor(1:prim_num) = prim_normalization_factor(1:prim_num)
-!    else
+    else
-!      factor(1:prim_num) = 1.d0
+      factor(1:prim_num) = 1.d0
-!    endif
+    endif
+
     rc = trexio_write_basis_prim_factor(f(1), factor)
     call trexio_assert(rc, TREXIO_SUCCESS)
     deallocate(factor)
@@ -324,14 +321,10 @@ subroutine export_trexio(update,full_path)
     C_A(3) = 0.d0
 
     allocate(factor(ao_num))
-    if (ao_normalized) then
+    do i=1,ao_num
-      do i=1,ao_num
+      l = ao_first_of_shell(ao_shell(i))
-        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))
-        factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
+    enddo
-      enddo
-    else
-      factor(:) = 1.d0
-    endif
     rc = trexio_write_ao_normalization(f(1), factor)
     call trexio_assert(rc, TREXIO_SUCCESS)
     deallocate(factor)