Fixes for numerical orbitals in qp_import

ocaml/Input_ao_basis.ml

@@ -44,8 +44,12 @@ end = struct
   let get_default = Qpackage.get_ezfio_default "ao_basis";;
 
   let read_ao_basis () =
-    Ezfio.get_ao_basis_ao_basis ()
-    |> AO_basis_name.of_string
+    let result =
+      Ezfio.get_ao_basis_ao_basis ()
+    in
+    if result <> "None" then
+      AO_basis_name.of_string result
+    else failwith "No basis"
   ;;
 
   let read_ao_num () =
@@ -192,7 +196,7 @@ end = struct
          ao_expo         ;
          ao_cartesian    ;
          ao_normalized   ;
-         primitives_normalized ; 
+         primitives_normalized ;
        } = b
      in
      write_md5 b ;
@@ -207,7 +211,7 @@ end = struct
      Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
        ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
 
-     let ao_nucl = 
+     let ao_nucl =
        Array.to_list ao_nucl
        |> list_map Nucl_number.to_int
      in
@@ -215,7 +219,7 @@ end = struct
        ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
 
      let ao_power =
-       let l = Array.to_list ao_power in 
+       let l = Array.to_list ao_power in
        List.concat [
          (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.x) l) ;
          (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.y) l) ;
@@ -227,7 +231,7 @@ end = struct
      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
       |> list_map AO_coef.to_float
@@ -267,7 +271,10 @@ end = struct
         |> Ezfio.set_ao_basis_ao_md5 ;
         Some result
       with
-      | _ -> (Ezfio.set_ao_basis_ao_md5 "None" ; None)
+      | _ -> ( "None"
+               |> Digest.string
+               |> Digest.to_hex
+               |> Ezfio.set_ao_basis_ao_md5 ; None)
   ;;
 
 
@@ -276,7 +283,7 @@ end = struct
       to_basis b
       |> Long_basis.of_basis
       |> Array.of_list
-    and unordered_basis = 
+    and unordered_basis =
       to_long_basis b
       |> Array.of_list
     in
@@ -289,15 +296,15 @@ end = struct
             (a.(i) <- None ; i)
           else
             find x a (i+1)
-      and find2 (s,g,n) a i = 
+      and find2 (s,g,n) a i =
         if i = Array.length a then -1
         else
-            match a.(i) with 
+            match a.(i) with
                 | None -> find2 (s,g,n) a (i+1)
                 | Some (s', g', n')  ->
                    if s <> s' || n <> n' then find2 (s,g,n) a (i+1)
                    else
-                   let lc  = list_map (fun (prim, _) -> prim) g.Gto.lc 
+                   let lc  = list_map (fun (prim, _) -> prim) g.Gto.lc
                    and lc' = list_map (fun (prim, _) -> prim) g'.Gto.lc
                    in
                    if lc <> lc' then find2 (s,g,n) a (i+1) else (a.(i) <- None ; i)
@@ -313,13 +320,13 @@ end = struct
       let ao_num = List.length long_basis |> AO_number.of_int in
       let ao_prim_num =
         list_map (fun (_,g,_) -> List.length g.Gto.lc
-          |> AO_prim_number.of_int ) long_basis 
+          |> AO_prim_number.of_int ) long_basis
         |> Array.of_list
       and ao_nucl =
-        list_map (fun (_,_,n) -> n) long_basis 
+        list_map (fun (_,_,n) -> n) long_basis
         |> Array.of_list
       and ao_power =
-        list_map (fun (x,_,_) -> x) long_basis 
+        list_map (fun (x,_,_) -> x) long_basis
         |> Array.of_list
       in
       let ao_prim_num_max = Array.fold_left (fun s x ->
@@ -329,16 +336,16 @@ end = struct
       in
 
       let gtos =
-        list_map (fun (_,x,_) -> x) long_basis 
+        list_map (fun (_,x,_) -> x) long_basis
       in
       let create_expo_coef ec =
           let coefs =
             begin match ec with
             | `Coefs -> list_map (fun x->
-              list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos 
+              list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos
             | `Expos -> list_map (fun x->
               list_map (fun (prim,_) -> AO_expo.to_float
-              prim.GaussianPrimitive.expo) x.Gto.lc ) gtos 
+              prim.GaussianPrimitive.expo) x.Gto.lc ) gtos
             end
           in
           let rec get_n n accu = function
@@ -360,7 +367,7 @@ end = struct
       let ao_coef = create_expo_coef `Coefs
       |> Array.of_list
       |> Array.map AO_coef.of_float
-      and ao_expo = create_expo_coef `Expos 
+      and ao_expo = create_expo_coef `Expos
       |> Array.of_list
       |> Array.map AO_expo.of_float
       in
@@ -372,7 +379,7 @@ end = struct
         }
   ;;
 
-  let reorder b = 
+  let reorder b =
     let order = ordering b in
     let f a = Array.init (Array.length a) (fun i -> a.(order.(i))) in
     let ao_prim_num_max = AO_prim_number.to_int b.ao_prim_num_max
@@ -464,7 +471,7 @@ Basis set (read-only) ::
     | line :: tail ->
       let line = String.trim line in
       if line = "Basis set (read-only) ::" then
-        String.concat "\n" tail 
+        String.concat "\n" tail
       else
         extract_basis tail
     in

ocaml/Input_mo_basis.ml

@@ -56,7 +56,10 @@ end = struct
   let read_ao_md5 () =
     let ao_md5 =
       match (Input_ao_basis.Ao_basis.read ()) with
-      | None -> failwith "Unable to read AO basis"
+      | None -> ("None"
+                 |> Digest.string
+                 |> Digest.to_hex
+                 |> MD5.of_string)
       | Some result -> Input_ao_basis.Ao_basis.to_md5 result
     in
     let result =

scripts/qp_import_trexio.py

@@ -132,60 +132,113 @@ def write_ezfio(trexio_filename, filename):
     try:
         basis_type = trexio.read_basis_type(trexio_file)
 
-        if basis_type.lower() not in ["gaussian", "slater"]:
-            raise TypeError
+        if basis_type.lower() in ["gaussian", "slater"]:
+            shell_num   = trexio.read_basis_shell_num(trexio_file)
+            prim_num    = trexio.read_basis_prim_num(trexio_file)
+            ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
+            nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
+            exponent    = trexio.read_basis_exponent(trexio_file)
+            coefficient = trexio.read_basis_coefficient(trexio_file)
+            shell_index = trexio.read_basis_shell_index(trexio_file)
+            ao_shell    = trexio.read_ao_shell(trexio_file)
 
-        shell_num   = trexio.read_basis_shell_num(trexio_file)
-        prim_num    = trexio.read_basis_prim_num(trexio_file)
-        ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
-        nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
-        exponent    = trexio.read_basis_exponent(trexio_file)
-        coefficient = trexio.read_basis_coefficient(trexio_file)
-        shell_index = trexio.read_basis_shell_index(trexio_file)
-        ao_shell    = trexio.read_ao_shell(trexio_file)
+            ezfio.set_basis_basis("Read from TREXIO")
+            ezfio.set_ao_basis_ao_basis("Read from TREXIO")
+            ezfio.set_basis_shell_num(shell_num)
+            ezfio.set_basis_prim_num(prim_num)
+            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)
 
-        ezfio.set_basis_basis("Read from TREXIO")
-        ezfio.set_basis_shell_num(shell_num)
-        ezfio.set_basis_prim_num(prim_num)
-        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
+            m = 0
+            for i in ao_shell:
+                if i != prev:
+                   m += 1
+                   if prev is None or nucl_index[i] != nucl_index[prev]:
+                        nucl_shell_num.append(m)
+                        m = 0
+                prev = i
+            assert (len(nucl_shell_num) == nucl_num)
 
-        nucl_shell_num = []
-        prev = None
-        m = 0
-        for i in ao_shell:
-            if i != prev:
-               m += 1
-               if prev is None or nucl_index[i] != nucl_index[prev]:
-                    nucl_shell_num.append(m)
-                    m = 0
-            prev = i
-        assert (len(nucl_shell_num) == nucl_num)
+            shell_prim_num = []
+            prev = shell_index[0]
+            count = 0
+            for i in shell_index:
+                if i != prev:
+                   shell_prim_num.append(count)
+                   count = 0
+                count += 1
+                prev = i
+            shell_prim_num.append(count)
 
-        shell_prim_num = []
-        prev = shell_index[0]
-        count = 0
-        for i in shell_index:
-            if i != prev:
-               shell_prim_num.append(count)
-               count = 0
-            count += 1
-            prev = i
-        shell_prim_num.append(count)
+            assert (len(shell_prim_num) == shell_num)
 
-        assert (len(shell_prim_num) == shell_num)
-
-        ezfio.set_basis_shell_prim_num(shell_prim_num)
-        ezfio.set_basis_shell_index([x+1 for x in shell_index])
-        ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
+            ezfio.set_basis_shell_prim_num(shell_prim_num)
+            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)
-        prim_factor  = trexio.read_basis_prim_factor(trexio_file)
+            shell_factor = trexio.read_basis_shell_factor(trexio_file)
+            prim_factor  = trexio.read_basis_prim_factor(trexio_file)
 
-        print("OK")
+        elif basis_type.lower() == "numerical":
+
+            shell_num   = trexio.read_basis_shell_num(trexio_file)
+            prim_num    = shell_num
+            ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
+            nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
+            exponent    = [1.]*prim_num
+            coefficient = [1.]*prim_num
+            shell_index = [i for i in range(shell_num)]
+            ao_shell    = trexio.read_ao_shell(trexio_file)
+
+            ezfio.set_basis_basis("None")
+            ezfio.set_ao_basis_ao_basis("None")
+            ezfio.set_basis_shell_num(shell_num)
+            ezfio.set_basis_prim_num(prim_num)
+            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
+            m = 0
+            for i in ao_shell:
+                if i != prev:
+                   m += 1
+                   if prev is None or nucl_index[i] != nucl_index[prev]:
+                        nucl_shell_num.append(m)
+                        m = 0
+                prev = i
+            assert (len(nucl_shell_num) == nucl_num)
+
+            shell_prim_num = []
+            prev = shell_index[0]
+            count = 0
+            for i in shell_index:
+                if i != prev:
+                   shell_prim_num.append(count)
+                   count = 0
+                count += 1
+                prev = i
+            shell_prim_num.append(count)
+
+            assert (len(shell_prim_num) == shell_num)
+
+            ezfio.set_basis_shell_prim_num(shell_prim_num)
+            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)
+            prim_factor  = [1.]*prim_num
+        else:
+           raise TypeError
+
+        print(basis_type)
     except:
         print("None")
         ezfio.set_ao_basis_ao_cartesian(True)
@@ -262,7 +315,6 @@ def write_ezfio(trexio_filename, filename):
 #        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)
-        ezfio.set_ao_basis_ao_basis("Read from TREXIO")
 
         print("OK")
 

src/trexio/import_trexio_integrals.irp.f

@@ -3,6 +3,7 @@ program import_integrals_ao
   implicit none
   integer(trexio_t)              :: f ! TREXIO file handle
   integer(trexio_exit_code)      :: rc
+  PROVIDE mo_num
 
   f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
   if (f == 0_8) then
@@ -42,10 +43,10 @@ subroutine run(f)
 
   if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
     rc = trexio_read_nucleus_repulsion(f, s)
-    call trexio_assert(rc, TREXIO_SUCCESS)
     if (rc /= TREXIO_SUCCESS) then
       print *, irp_here, rc
       print *, 'Error reading nuclear repulsion'
+      call trexio_assert(rc, TREXIO_SUCCESS)
       stop -1
     endif
     call ezfio_set_nuclei_nuclear_repulsion(s)
@@ -63,6 +64,7 @@ subroutine run(f)
     if (rc /= TREXIO_SUCCESS) then
       print *, irp_here
       print *, 'Error reading AO overlap'
+      call trexio_assert(rc, TREXIO_SUCCESS)
       stop -1
     endif
     call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A)
@@ -74,6 +76,7 @@ subroutine run(f)
     if (rc /= TREXIO_SUCCESS) then
       print *, irp_here
       print *, 'Error reading AO kinetic integrals'
+      call trexio_assert(rc, TREXIO_SUCCESS)
       stop -1
     endif
     call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A)
@@ -85,6 +88,7 @@ subroutine run(f)
 !    if (rc /= TREXIO_SUCCESS) then
 !      print *, irp_here
 !      print *, 'Error reading AO ECP local integrals'
+!      call trexio_assert(rc, TREXIO_SUCCESS)
 !      stop -1
 !    endif
 !    call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A)
@@ -96,6 +100,7 @@ subroutine run(f)
     if (rc /= TREXIO_SUCCESS) then
       print *, irp_here
       print *, 'Error reading AO potential N-e integrals'
+      call trexio_assert(rc, TREXIO_SUCCESS)
       stop -1
     endif
     call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A)
@@ -106,41 +111,112 @@ subroutine run(f)
 
   ! AO 2e integrals
   ! ---------------
-  PROVIDE ao_integrals_map
 
-  integer*4 :: BUFSIZE
-  BUFSIZE=ao_num**2
-  allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
-  allocate(Vi(4,BUFSIZE), V(BUFSIZE))
+  rc = trexio_has_ao_2e_int(f)
+  PROVIDE ao_num
+  if (rc /= TREXIO_HAS_NOT) then
+      PROVIDE ao_integrals_map
 
-  integer*8 :: offset, icount
+      integer*4 :: BUFSIZE
+      BUFSIZE=ao_num**2
+      allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+      allocate(Vi(4,BUFSIZE), V(BUFSIZE))
 
-  offset = 0_8
-  icount = BUFSIZE
-  rc = TREXIO_SUCCESS
-  do while (icount == size(V))
-    rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
-    do m=1,icount
-      i = Vi(1,m)
-      j = Vi(2,m)
-      k = Vi(3,m)
-      l = Vi(4,m)
-      integral = V(m)
-      call two_e_integrals_index(i, j, k, l, buffer_i(m) )
-      buffer_values(m) = integral
-    enddo
-    call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
-    offset = offset + icount
+      integer*8 :: offset, icount
+
+      offset = 0_8
+      icount = BUFSIZE
+      rc = TREXIO_SUCCESS
+      do while (icount == size(V))
+        rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
+        do m=1,icount
+          i = Vi(1,m)
+          j = Vi(2,m)
+          k = Vi(3,m)
+          l = Vi(4,m)
+          integral = V(m)
+          call two_e_integrals_index(i, j, k, l, buffer_i(m) )
+          buffer_values(m) = integral
+        enddo
+        call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
+        offset = offset + icount
+        if (rc /= TREXIO_SUCCESS) then
+            exit
+        endif
+      end do
+      n_integrals = offset
+
+      call map_sort(ao_integrals_map)
+      call map_unique(ao_integrals_map)
+
+      call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+      call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+
+      deallocate(buffer_i, buffer_values, Vi, V)
+      print *, 'AO integrals read from TREXIO file'
+  else
+      print *, 'AO integrals not found in TREXIO file'
+  endif
+
+  ! MO integrals
+  ! ------------
+
+  allocate(A(mo_num, mo_num))
+  if (trexio_has_mo_1e_int_core_hamiltonian(f) == TREXIO_SUCCESS) then
+    rc = trexio_read_mo_1e_int_core_hamiltonian(f, A)
     if (rc /= TREXIO_SUCCESS) then
-        exit
+      print *, irp_here
+      print *, 'Error reading MO 1e integrals'
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      stop -1
     endif
-  end do
-  n_integrals = offset
+    call ezfio_set_mo_one_e_ints_mo_one_e_integrals(A)
+    call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('Read')
+  endif
+  deallocate(A)
 
-  call map_sort(ao_integrals_map)
-  call map_unique(ao_integrals_map)
+  ! MO 2e integrals
+  ! ---------------
 
-  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
-  call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+  rc = trexio_has_mo_2e_int(f)
+  if (rc /= TREXIO_HAS_NOT) then
+
+      BUFSIZE=mo_num**2
+      allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+      allocate(Vi(4,BUFSIZE), V(BUFSIZE))
+
+
+      offset = 0_8
+      icount = BUFSIZE
+      rc = TREXIO_SUCCESS
+      do while (icount == size(V))
+        rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V)
+        do m=1,icount
+          i = Vi(1,m)
+          j = Vi(2,m)
+          k = Vi(3,m)
+          l = Vi(4,m)
+          integral = V(m)
+          call two_e_integrals_index(i, j, k, l, buffer_i(m) )
+          buffer_values(m) = integral
+        enddo
+        call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4))
+        offset = offset + icount
+        if (rc /= TREXIO_SUCCESS) then
+            exit
+        endif
+      end do
+      n_integrals = offset
+
+      call map_sort(mo_integrals_map)
+      call map_unique(mo_integrals_map)
+
+      call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
+      call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
+      deallocate(buffer_i, buffer_values, Vi, V)
+      print *, 'MO integrals read from TREXIO file'
+  else
+      print *, 'MO integrals not found in TREXIO file'
+  endif
 
 end