Improvements to read cholesky integrals from a trexio file

ocaml/Input_ao_basis.ml

@@ -58,52 +58,70 @@ end = struct
   ;;
 
   let read_ao_prim_num () =
-    Ezfio.get_ao_basis_ao_prim_num ()
+    if Ezfio.has_ao_basis_ao_prim_num () then
-    |> Ezfio.flattened_ezfio
+      Ezfio.get_ao_basis_ao_prim_num ()
-    |> Array.map AO_prim_number.of_int
+      |> Ezfio.flattened_ezfio
+      |> Array.map AO_prim_number.of_int
+    else
+      [||]
   ;;
 
   let read_ao_prim_num_max () =
-    Ezfio.get_ao_basis_ao_prim_num ()
+    if Ezfio.has_ao_basis_ao_prim_num () then
-    |> Ezfio.flattened_ezfio
+      Ezfio.get_ao_basis_ao_prim_num ()
-    |> Array.fold_left (fun x y -> if x>y then x else y) 0
+      |> Ezfio.flattened_ezfio
-    |> AO_prim_number.of_int
+      |> Array.fold_left (fun x y -> if x>y then x else y) 0
+      |> AO_prim_number.of_int
+    else
+      AO_prim_number.of_int 0
   ;;
 
   let read_ao_nucl () =
-    let nmax = Nucl_number.get_max () in
+    if Ezfio.has_ao_basis_ao_nucl () then
-    Ezfio.get_ao_basis_ao_nucl ()
+      let nmax = Nucl_number.get_max () in
-    |> Ezfio.flattened_ezfio
+      Ezfio.get_ao_basis_ao_nucl ()
-    |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x)
+      |> Ezfio.flattened_ezfio
+      |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x)
+    else
+      [||]
   ;;
 
   let read_ao_power () =
-    let x = Ezfio.get_ao_basis_ao_power () in
+    if Ezfio.has_ao_basis_ao_power () then
-    let dim = x.Ezfio.dim.(0) in
+      let x = Ezfio.get_ao_basis_ao_power () in
-    let data = Ezfio.flattened_ezfio x in
+      let dim = x.Ezfio.dim.(0) in
-    let result = Array.init dim (fun x -> "") in
+      let data = Ezfio.flattened_ezfio x in
-    for i=1 to dim
+      let result = Array.init dim (fun x -> "") in
-    do
+      for i=1 to dim
-      if (data.(i-1) > 0) then
+      do
-        result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1));
+        if (data.(i-1) > 0) then
-      if (data.(dim+i-1) > 0) then
+          result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1));
-        result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1));
+        if (data.(dim+i-1) > 0) then
-      if (data.(2*dim+i-1) > 0) then
+          result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1));
-        result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1));
+        if (data.(2*dim+i-1) > 0) then
-    done;
+          result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1));
-    Array.map Angmom.Xyz.of_string result
+      done;
+      Array.map Angmom.Xyz.of_string result
+    else
+      [||]
   ;;
 
   let read_ao_coef () =
-    Ezfio.get_ao_basis_ao_coef ()
+    if Ezfio.has_ao_basis_ao_coef () then
-    |> Ezfio.flattened_ezfio
+      Ezfio.get_ao_basis_ao_coef ()
-    |> Array.map AO_coef.of_float
+      |> Ezfio.flattened_ezfio
+      |> Array.map AO_coef.of_float
+    else
+      [||]
   ;;
 
   let read_ao_expo () =
-    Ezfio.get_ao_basis_ao_expo ()
+    if Ezfio.has_ao_basis_ao_expo () then
-    |> Ezfio.flattened_ezfio
+      Ezfio.get_ao_basis_ao_expo ()
-    |> Array.map AO_expo.of_float
+      |> Ezfio.flattened_ezfio
+      |> Array.map AO_expo.of_float
+    else
+      [||]
   ;;
 
   let read_ao_cartesian () =

src/ao_two_e_ints/cholesky.irp.f

@@ -79,7 +79,7 @@ END_PROVIDER
    type(mmap_type) :: map
 
    PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem
-   PROVIDE nucl_coord ao_two_e_integral_schwartz
+   PROVIDE nucl_coord
    call set_multiple_levels_omp(.False.)
 
    call wall_time(wall0)

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -1,6 +1,22 @@
 
 ! ---
 
+logical function do_schwartz_accel(i,j,k,l)
+ implicit none
+ BEGIN_DOC
+ ! If true, use Schwatrz to accelerate direct integral calculation
+ END_DOC
+  integer, intent(in)            :: i, j, k, l
+  if (do_ao_cholesky) then
+    do_schwartz_accel = .False.
+  else
+    do_schwartz_accel = (ao_prim_num(i) * ao_prim_num(j) * &
+                         ao_prim_num(k) * ao_prim_num(l) > 1024 )
+  endif
+
+end function
+
+
 double precision function ao_two_e_integral(i, j, k, l)
 
   BEGIN_DOC
@@ -25,6 +41,7 @@ double precision function ao_two_e_integral(i, j, k, l)
   double precision, external     :: ao_two_e_integral_cosgtos
   double precision, external     :: ao_two_e_integral_schwartz_accel
 
+  logical, external :: do_schwartz_accel
 
   if(use_cosgtos) then
     !print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
@@ -35,7 +52,7 @@ double precision function ao_two_e_integral(i, j, k, l)
 
     ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l)
 
-  else if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
+  else if (do_schwartz_accel(i,j,k,l)) then
 
        ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
 

src/tools/print_energy.irp.f

@@ -16,20 +16,20 @@ subroutine run
  implicit none
  call print_mol_properties
  print *, psi_energy + nuclear_repulsion
-! call print_energy_components
+ call print_energy_components
-! print *, 'E(HF) = ', HF_energy
+ print *, 'E(HF) = ', HF_energy
-! print *, 'E(CI) = ', psi_energy + nuclear_repulsion
+ print *, 'E(CI) = ', psi_energy + nuclear_repulsion
-! print *, ''
+ print *, ''
-! print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy
+ print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy
-! print *, 'E_kin(HF) = ', HF_kinetic_energy
+ print *, 'E_kin(HF) = ', HF_kinetic_energy
-! print *, ''
+ print *, ''
-! print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy
+ print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy
-! print *, 'E_ne (HF) = ', HF_n_e_energy
+ print *, 'E_ne (HF) = ', HF_n_e_energy
-! print *, ''
+ print *, ''
-! print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy
+ print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy
-! print *, 'E_1e (HF) = ', HF_one_electron_energy
+ print *, 'E_1e (HF) = ', HF_one_electron_energy
-! print *, ''
+ print *, ''
-! print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy
+ print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy
-! print *, 'E_2e (HF) = ', HF_two_electron_energy
+ print *, 'E_2e (HF) = ', HF_two_electron_energy
 
 end

src/trexio/import_trexio_integrals.irp.f

@@ -46,6 +46,8 @@ subroutine run(f)
   double precision, allocatable :: tmp(:,:,:)
   integer*8 :: offset, icount
 
+  integer :: k_num
+
   integer, external :: getUnitAndOpen
 
   if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
@@ -163,46 +165,48 @@ subroutine run(f)
 
           deallocate(Vi, V, tmp)
           print *, 'Cholesky AO integrals read from TREXIO file'
-      endif
 
-      rc = trexio_has_ao_2e_int_eri(f)
+      else
-      if (rc /= TREXIO_HAS_NOT) then
-          PROVIDE ao_integrals_map
 
-          BUFSIZE=ao_num**2
+          rc = trexio_has_ao_2e_int_eri(f)
-          allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+          if (rc /= TREXIO_HAS_NOT) then
-          allocate(Vi(4,BUFSIZE), V(BUFSIZE))
+              PROVIDE ao_integrals_map
 
-          offset = 0_8
+              BUFSIZE=ao_num**2
-          icount = BUFSIZE
+              allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
-          rc = TREXIO_SUCCESS
+              allocate(Vi(4,BUFSIZE), V(BUFSIZE))
-          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)
+              offset = 0_8
-          call map_unique(ao_integrals_map)
+              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_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+              call map_sort(ao_integrals_map)
-          call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+              call map_unique(ao_integrals_map)
 
-          deallocate(buffer_i, buffer_values, Vi, V)
+              call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
-          print *, 'AO integrals read from TREXIO file'
+              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'
+          endif
       endif
   else
       print *, 'AO integrals not found in TREXIO file'
@@ -270,44 +274,47 @@ subroutine run(f)
 
           deallocate(Vi, V, tmp)
           print *, 'Cholesky MO integrals read from TREXIO file'
-      endif
 
-      rc = trexio_has_mo_2e_int_eri(f)
+      else
-      if (rc /= TREXIO_HAS_NOT) then
+
-          BUFSIZE=mo_num**2
+        rc = trexio_has_mo_2e_int_eri(f)
-          allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+        if (rc /= TREXIO_HAS_NOT) then
-          allocate(Vi(4,BUFSIZE), V(BUFSIZE))
+            BUFSIZE=mo_num**2
+            allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+            allocate(Vi(4,BUFSIZE), V(BUFSIZE))
 
 
-          offset = 0_8
+            offset = 0_8
-          icount = BUFSIZE
+            icount = BUFSIZE
-          rc = TREXIO_SUCCESS
+            rc = TREXIO_SUCCESS
-          do while (icount == size(V))
+            do while (icount == size(V))
-            rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V)
+              rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V)
-            do m=1,icount
+              do m=1,icount
-              i = Vi(1,m)
+                i = Vi(1,m)
-              j = Vi(2,m)
+                j = Vi(2,m)
-              k = Vi(3,m)
+                k = Vi(3,m)
-              l = Vi(4,m)
+                l = Vi(4,m)
-              integral = V(m)
+                integral = V(m)
-              call two_e_integrals_index(i, j, k, l, buffer_i(m) )
+                call two_e_integrals_index(i, j, k, l, buffer_i(m) )
-              buffer_values(m) = integral
+                buffer_values(m) = integral
-            enddo
+              enddo
-            call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4))
+              call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4))
-            offset = offset + icount
+              offset = offset + icount
-            if (rc /= TREXIO_SUCCESS) then
+              if (rc /= TREXIO_SUCCESS) then
-                exit
+                  exit
-            endif
+              endif
-          end do
+            end do
-          n_integrals = offset
+            n_integrals = offset
 
-          call map_sort(mo_integrals_map)
+            call map_sort(mo_integrals_map)
-          call map_unique(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'
+        endif
 
-          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'
       endif
 
   else