From b3930d977ac6d6cfd1bfb1dcf76a1b23e1e8b090 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 4 Oct 2024 21:52:05 +0200 Subject: [PATCH] Improvements to read cholesky integrals from a trexio file --- ocaml/Input_ao_basis.ml | 80 ++++++++----- src/ao_two_e_ints/cholesky.irp.f | 2 +- src/ao_two_e_ints/two_e_integrals.irp.f | 19 ++- src/tools/print_energy.irp.f | 30 ++--- src/trexio/import_trexio_integrals.irp.f | 141 ++++++++++++----------- 5 files changed, 157 insertions(+), 115 deletions(-) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index 506cf069..d9e28e04 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -58,52 +58,70 @@ end = struct ;; let read_ao_prim_num () = - Ezfio.get_ao_basis_ao_prim_num () - |> Ezfio.flattened_ezfio - |> Array.map AO_prim_number.of_int + if Ezfio.has_ao_basis_ao_prim_num () then + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> Array.map AO_prim_number.of_int + else + [||] ;; let read_ao_prim_num_max () = - Ezfio.get_ao_basis_ao_prim_num () - |> Ezfio.flattened_ezfio - |> Array.fold_left (fun x y -> if x>y then x else y) 0 - |> AO_prim_number.of_int + if Ezfio.has_ao_basis_ao_prim_num () then + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> 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 - Ezfio.get_ao_basis_ao_nucl () - |> Ezfio.flattened_ezfio - |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x) + if Ezfio.has_ao_basis_ao_nucl () then + let nmax = Nucl_number.get_max () in + Ezfio.get_ao_basis_ao_nucl () + |> 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 - let dim = x.Ezfio.dim.(0) in - let data = Ezfio.flattened_ezfio x in - let result = Array.init dim (fun x -> "") in - for i=1 to dim - do - if (data.(i-1) > 0) then - result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1)); - if (data.(dim+i-1) > 0) then - result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1)); - if (data.(2*dim+i-1) > 0) then - result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); - done; - Array.map Angmom.Xyz.of_string result + if Ezfio.has_ao_basis_ao_power () then + let x = Ezfio.get_ao_basis_ao_power () in + let dim = x.Ezfio.dim.(0) in + let data = Ezfio.flattened_ezfio x in + let result = Array.init dim (fun x -> "") in + for i=1 to dim + do + if (data.(i-1) > 0) then + result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1)); + if (data.(dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1)); + if (data.(2*dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); + done; + Array.map Angmom.Xyz.of_string result + else + [||] ;; let read_ao_coef () = - Ezfio.get_ao_basis_ao_coef () - |> Ezfio.flattened_ezfio - |> Array.map AO_coef.of_float + if Ezfio.has_ao_basis_ao_coef () then + Ezfio.get_ao_basis_ao_coef () + |> Ezfio.flattened_ezfio + |> Array.map AO_coef.of_float + else + [||] ;; let read_ao_expo () = - Ezfio.get_ao_basis_ao_expo () - |> Ezfio.flattened_ezfio - |> Array.map AO_expo.of_float + if Ezfio.has_ao_basis_ao_expo () then + Ezfio.get_ao_basis_ao_expo () + |> Ezfio.flattened_ezfio + |> Array.map AO_expo.of_float + else + [||] ;; let read_ao_cartesian () = diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index d409a6d5..efafd504 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/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) diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index d12f3d45..7b4a2e5a 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/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) diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f index 1d3ea648..485a0910 100644 --- a/src/tools/print_energy.irp.f +++ b/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 -! print *, 'E(HF) = ', HF_energy -! print *, 'E(CI) = ', psi_energy + nuclear_repulsion -! print *, '' -! print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy -! print *, 'E_kin(HF) = ', HF_kinetic_energy -! print *, '' -! print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy -! print *, 'E_ne (HF) = ', HF_n_e_energy -! print *, '' -! print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy -! print *, 'E_1e (HF) = ', HF_one_electron_energy -! print *, '' -! print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy -! print *, 'E_2e (HF) = ', HF_two_electron_energy + call print_energy_components + print *, 'E(HF) = ', HF_energy + print *, 'E(CI) = ', psi_energy + nuclear_repulsion + print *, '' + print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy + print *, 'E_kin(HF) = ', HF_kinetic_energy + print *, '' + print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy + print *, 'E_ne (HF) = ', HF_n_e_energy + print *, '' + print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy + print *, 'E_1e (HF) = ', HF_one_electron_energy + print *, '' + print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy + print *, 'E_2e (HF) = ', HF_two_electron_energy end diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f index 556ed7bc..405c0863 100644 --- a/src/trexio/import_trexio_integrals.irp.f +++ b/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) - if (rc /= TREXIO_HAS_NOT) then - PROVIDE ao_integrals_map + else - BUFSIZE=ao_num**2 - allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) - allocate(Vi(4,BUFSIZE), V(BUFSIZE)) + rc = trexio_has_ao_2e_int_eri(f) + if (rc /= TREXIO_HAS_NOT) then + PROVIDE ao_integrals_map - 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 + BUFSIZE=ao_num**2 + allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) + allocate(Vi(4,BUFSIZE), V(BUFSIZE)) - call map_sort(ao_integrals_map) - call map_unique(ao_integrals_map) + 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_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') + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) - deallocate(buffer_i, buffer_values, Vi, V) - print *, 'AO integrals read from TREXIO file' + 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' + 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) - if (rc /= TREXIO_HAS_NOT) then - BUFSIZE=mo_num**2 - allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) - allocate(Vi(4,BUFSIZE), V(BUFSIZE)) + else + + rc = trexio_has_mo_2e_int_eri(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 + 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_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' + 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