From 0b29bb97d1ba19b53957ba11f53cbbe3b07dbbbb Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert Date: Thu, 24 Apr 2025 11:14:22 +0200 Subject: [PATCH] Use npz format to store post-processing results --- app/cli.f90 | 6 + app/cli_help.f90 | 3 + app/driver_run.f90 | 5 + man/tblite-run.1.adoc | 11 + src/tblite/double_dictionary.f90 | 630 ++++++++++++++------------- src/tblite/io/numpy/CMakeLists.txt | 1 + src/tblite/io/numpy/crc32.f90 | 4 +- src/tblite/io/numpy/meson.build | 1 + src/tblite/io/numpy/savez.f90 | 217 ++------- src/tblite/io/numpy/zip.f90 | 191 ++++++++ src/tblite/xtb/singlepoint.f90 | 1 - test/unit/test_double_dictionary.f90 | 145 +++--- test/unit/test_xtbml.f90 | 61 ++- 13 files changed, 678 insertions(+), 598 deletions(-) create mode 100644 src/tblite/io/numpy/zip.f90 diff --git a/app/cli.f90 b/app/cli.f90 index 49335ca6..f2a51244 100644 --- a/app/cli.f90 +++ b/app/cli.f90 @@ -72,6 +72,8 @@ module tblite_cli type(solvation_input), allocatable :: solvation !> Input for post processing container character(len=:), allocatable :: post_processing + !> Input for post processing container + character(len=:), allocatable :: post_proc_output !> Numerical accuracy for self-consistent iterations real(wp) :: accuracy = 1.0_wp !> Maximum number of iterations for SCF @@ -567,6 +569,10 @@ subroutine get_run_arguments(config, list, start, error) end if end if + if (allocated(config%post_processing) .and. .not.allocated(config%post_proc_output)) then + config%post_proc_output = "tblite-data.npz" + end if + if (allocated(solvent)) then if (.not.allocated(sol_state)) then sol_state = solution_state%gsolv diff --git a/app/cli_help.f90 b/app/cli_help.f90 index 437a2163..0697b9a6 100644 --- a/app/cli_help.f90 +++ b/app/cli_help.f90 @@ -150,6 +150,9 @@ module tblite_cli_help " Mayer-Wiberg bond orders are computed by default."//nl//& " Options: molmom (molecular multipole moments)"//nl//& " Options: xtbml (atomistic properties based on Mulliken partitioning)"//nl//& + " --post-processing-output "//nl//& + " Output file for post processing results in npz format."//nl//& + " (default: tblite-data.npz)"//nl//& " --grad [file] Evaluate molecular gradient and virial"//nl//& " Results are stored in file (default: tblite.txt)"//nl//& " --json [file] Dump results as JSON output (default: tblite.json)"//nl//& diff --git a/app/driver_run.f90 b/app/driver_run.f90 index 2abb2f6e..73f7f948 100644 --- a/app/driver_run.f90 +++ b/app/driver_run.f90 @@ -302,6 +302,11 @@ subroutine run_main(config, error) if (allocated(error)) return end if + if (allocated(post_proc) .and. allocated(config%post_proc_output)) then + call results%dict%dump(config%post_proc_output, error) + if (allocated(error)) return + end if + if (config%verbosity > 2) then call ascii_levels(ctx%unit, config%verbosity, wfn%homo, wfn%emo, wfn%focc, 7) call post_proc%dict%get_entry("molecular-dipole", dpmom) diff --git a/man/tblite-run.1.adoc b/man/tblite-run.1.adoc index 75c4e5d6..ab04bbf5 100644 --- a/man/tblite-run.1.adoc +++ b/man/tblite-run.1.adoc @@ -112,6 +112,17 @@ Supported geometry input formats are: Evaluates analytical gradient, results are stored in _file_ (default: tblite.txt) +*--post-processing* _file_|_name_:: + Add post processing methods to the calculation either by using a toml file as input + or by specifying the post-processing step name directly. + Mayer-Wiberg bond orders are computed by default. + - molmom (molecular multipole moments) + - xtbml (atomistic properties based on Mulliken partitioning) + +*--post-processing-output* _file_:: + Output file for post processing results in npz format. + (default: tblite-data.npz) + *--json* [_file_]:: Dump results as JSON output to _file_ (default: tblite.json) diff --git a/src/tblite/double_dictionary.f90 b/src/tblite/double_dictionary.f90 index 2dfd9f89..1194af73 100644 --- a/src/tblite/double_dictionary.f90 +++ b/src/tblite/double_dictionary.f90 @@ -24,6 +24,8 @@ module tblite_double_dictionary use mctc_env, only : error_type, fatal_error use tblite_toml, only : toml_array, toml_table, toml_key, add_table, set_value, toml_error use tblite_toml, only : toml_dump, add_array, get_value, toml_parse + use tblite_io_numpy, only : save_npz, load_npz + use tblite_io_numpy_zip, only : zip_file, list_zip_file implicit none private @@ -81,23 +83,24 @@ module tblite_double_dictionary procedure :: remove_entry_index generic, public :: get_index => return_label_index procedure :: return_label_index - procedure :: dump_to_toml procedure :: dump_to_file - procedure :: dump_to_unit - generic, public :: dump => dump_to_file, dump_to_toml, dump_to_unit + generic, public :: dump => dump_to_file generic, public :: operator(==) => equal_dict procedure :: equal_dict - generic, public :: load => load_from_file, load_from_unit, load_from_toml - procedure :: load_from_toml + generic, public :: load => load_from_file procedure :: load_from_file - procedure :: load_from_unit end type double_dictionary_type contains -function equal_record(lhs, rhs) result(equal) - class(double_record), intent(in) :: lhs, rhs +!> Check if two double records are equal +pure function equal_record(lhs, rhs) result(equal) + !> First record to compare + class(double_record), intent(in) :: lhs + !> Second record to compare + class(double_record), intent(in) :: rhs + !> Result of the comparison logical :: equal equal = .false. @@ -106,178 +109,116 @@ function equal_record(lhs, rhs) result(equal) end if if (allocated(lhs%array1) .and. allocated(rhs%array1)) then - if (all(lhs%array1 == rhs%array1)) then - equal = .true. - return - else - return - end if - endif + equal = all(lhs%array1 == rhs%array1) + return + end if if (allocated(lhs%array2) .and. allocated(rhs%array2)) then - if (all(lhs%array2 == rhs%array2)) then - equal = .true. - return - else - return - end if - endif + equal = all(lhs%array2 == rhs%array2) + return + end if if (allocated(lhs%array3) .and. allocated(rhs%array3)) then - if (all(lhs%array3 == rhs%array3)) then - equal = .true. - return - else - return - end if - endif + equal = all(lhs%array3 == rhs%array3) + return + end if -end function +end function equal_record -function equal_dict(lhs, rhs) result(equal) - class(double_dictionary_type), intent(in) :: lhs, rhs - integer :: i +!> Check if two double dictionaries are equal +pure function equal_dict(lhs, rhs) result(equal) + !> First dictionary to compare + class(double_dictionary_type), intent(in) :: lhs + !> Second dictionary to compare + class(double_dictionary_type), intent(in) :: rhs + !> Result of the comparison logical :: equal + + integer :: i equal = .false. - if (lhs%get_n_entries() /= rhs%get_n_entries()) then + if (lhs%n /= rhs%n) then return end if - do i = 1, lhs%get_n_entries() + do i = 1, lhs%n if (.not.(lhs%record(i) == rhs%record(i))) return end do equal = .true. -end function +end function equal_dict !> Read double dictionary data from file -subroutine load_from_file(self, file, error) +subroutine load_from_file(self, filename, error) !> Instance of the parametrization data class(double_dictionary_type), intent(inout) :: self !> File name - character(len=*), intent(in) :: file + character(len=*), intent(in) :: filename !> Error handling type(error_type), allocatable, intent(out) :: error - integer :: unit + integer :: io, stat, irec, it logical :: exist + type(zip_file) :: zip + character(len=:), allocatable :: msg, label + character(len=*), parameter :: prefix = "tblite0_" - inquire(file=file, exist=exist) - if (.not.exist) then - call fatal_error(error, "Could not find toml file '"//file//"'") - return + inquire(file=filename, exist=exist) + if (.not. exist) then + call fatal_error(error, "File '"//filename//"' does not exist") + return end if - open(file=file, newunit=unit) - call self%load(unit, error) - close(unit) -end subroutine load_from_file - - -!> Read double_dictionary data from file -subroutine load_from_unit(self, unit, error) - !> Instance of the double dictionary data - class(double_dictionary_type), intent(inout) :: self - !> File name - integer, intent(in) :: unit - !> Error handling - type(error_type), allocatable, intent(out) :: error - - type(toml_error), allocatable :: parse_error - type(toml_table), allocatable :: table - - call toml_parse(table, unit, parse_error) - - if (allocated(parse_error)) then - allocate(error) - call move_alloc(parse_error%message, error%message) + open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) + call list_zip_file(io, filename, zip, stat, msg) + close(io) + if (stat /= 0) then + if (.not.allocated(msg)) msg = "Failed to read zip file '"//filename//"'" + call fatal_error(error, msg) return end if - call self%load(table, error) - if (allocated(error)) return - -end subroutine load_from_unit - -subroutine load_from_toml(self, table, error) - use tblite_toml, only : len - !iterate over entries and dump to toml - class(double_dictionary_type) :: self - !> toml table to add entries to - type(toml_table), intent(inout) :: table - type(toml_key), allocatable :: list_keys(:) - !> Error handling - type(error_type), allocatable, intent(out) :: error - - type(toml_array), pointer :: array - real(kind=wp), allocatable :: array1(:) - - integer :: i, stat - - call table%get_keys(list_keys) - - do i = 1, size(list_keys) - call get_value(table, list_keys(i), array, stat=stat) - if (stat /= 0) then - call fatal_error(error, "Cannot read entry for array") - return - end if - if (allocated(array1)) deallocate(array1) - allocate(array1(len(array))) - call get_value(array, array1) - call self%add_entry(list_keys(i)%key, array1) - end do -end subroutine - -subroutine dump_to_toml(self, table, error) - !iterate over entries and dump to toml - class(double_dictionary_type) :: self - !> toml table to add entries to - type(toml_table), intent(inout) :: table - !> Error handling - type(error_type), allocatable, intent(out) :: error - type(toml_array), pointer :: array - real(kind=wp), allocatable :: array1(:), array2(:, :), array3(:, :, :) - - integer :: i, stat, j - do i = 1, self%get_n_entries() - call add_array(table, self%record(i)%label, array) - - if (allocated(array1)) deallocate(array1) - call self%get_entry(i, array1) - if (allocated(array1)) then - do j = 1, size(array1) - call set_value(array, j, array1(j), stat=stat) - end do - cycle - end if - - call self%get_entry(i, array2) - if (allocated(array2)) then - array1 = reshape(array2, [size(array2, 1)*size(array2, 2)]) - deallocate(array2) - do j = 1, size(array1) - call set_value(array, j, array1(j), stat=stat) - end do - cycle - end if + if (.not.allocated(zip%records)) then + call fatal_error(error, "No records found in file '"//filename//"'") + return + end if - call self%get_entry(i, array3) - if (allocated(array3)) then - array1 = reshape(array3, [size(array3, 1)*size(array3, 2)*size(array3, 3)]) - do j = 1, size(array1) - call set_value(array, j, array1(j), stat=stat) - end do - deallocate(array3) - cycle - end if + do irec = 1, size(zip%records) + associate(path => zip%records(irec)%path) + if (len(path) < len(prefix) + 3 + 4) then + call fatal_error(error, "File '"//filename//"::"//path// & + & "' does not match expected format") + exit + end if + label = path(len(prefix) + 3 + 1:len(path) - 4) + select case(path(:len(prefix) + 3)) + case default + call fatal_error(error, "Unknown file type '"//filename//"::"//path//"'") + exit + case (prefix//"r1_") + call self%push(label, it) + call load_npz(filename, prefix//"r1_"//label, self%record(it)%array1, stat, msg) + case (prefix//"r2_") + call self%push(label, it) + call load_npz(filename, prefix//"r2_"//label, self%record(it)%array2, stat, msg) + case (prefix//"r3_") + call self%push(label, it) + call load_npz(filename, prefix//"r3_"//label, self%record(it)%array3, stat, msg) + end select + if (stat /= 0) then + if (.not.allocated(msg)) then + msg = "Failed to load file '"//filename//"::"//path//"'" + end if + call fatal_error(error, msg) + exit + end if + end associate end do - deallocate(array1) -end subroutine +end subroutine load_from_file +!> Write double dictionary data to file subroutine dump_to_file(self, file, error) !> Instance of the parametrization data class(double_dictionary_type), intent(in) :: self @@ -286,87 +227,92 @@ subroutine dump_to_file(self, file, error) !> Error handling type(error_type), allocatable, intent(out) :: error - integer :: unit + character(len=*), parameter :: prefix = "tblite0_" - open(file=file, newunit=unit) - call self%dump(unit, error) - close(unit) - if (allocated(error)) return - -end subroutine dump_to_file + integer :: it, stat, io + logical :: exist + real(kind=wp), allocatable :: array1(:), array2(:, :), array3(:, :, :) + character(len=:), allocatable :: msg -!> Write double dictionary data to file -subroutine dump_to_unit(self, unit, error) - !> Instance of the parametrization data - class(double_dictionary_type), intent(in) :: self - !> Formatted unit - integer, intent(in) :: unit - !> Error handling - type(error_type), allocatable, intent(out) :: error + inquire(file=file, exist=exist) + if (exist) then + open(file=file, newunit=io) + close(io, status='delete') + end if - type(toml_table) :: table - type(toml_error), allocatable :: ser_error + do it = 1, self%n + associate(record => self%record(it)) + if (allocated(record%array1)) then + call save_npz(file, prefix//"r1_"//record%label, record%array1, stat, msg) + if (stat /= 0) exit + cycle + end if - table = toml_table() + if (allocated(record%array2)) then + call save_npz(file, prefix//"r2_"//record%label, record%array2, stat, msg) + if (stat /= 0) exit + cycle + end if - call self%dump(table, error) + if (allocated(record%array3)) then + call save_npz(file, prefix//"r3_"//record%label, record%array3, stat, msg) + if (stat /= 0) exit + cycle + end if + end associate + end do - call toml_dump(table, unit, ser_error) - if (allocated(ser_error)) then - call fatal_error(error, ser_error%message) - end if + if (allocated(error)) return +end subroutine dump_to_file -end subroutine dump_to_unit -subroutine remove_entry_index(self, index) - class(double_dictionary_type) :: self - integer :: index, old_n, i, it - type(double_dictionary_type) :: tmp +pure subroutine remove_entry_index(self, index) + class(double_dictionary_type), intent(inout) :: self + integer, intent(in) :: index + integer :: old_n, i, it + type(double_record), allocatable :: tmp(:) if (index > self%n) return - tmp = self + call move_alloc(self%record, tmp) old_n = self%n self%n = self%n - 1 - deallocate(self%record) allocate(self%record(self%n)) it = 1 do i = 1, old_n if (i == index) cycle - self%record(it) = tmp%record(i) + call move_record(tmp(i), self%record(it)) it = it + 1 end do -end subroutine +end subroutine remove_entry_index -subroutine remove_entry_label(self, label) - class(double_dictionary_type) :: self - character(len=*) :: label +pure subroutine remove_entry_label(self, label) + class(double_dictionary_type), intent(inout) :: self + character(len=*), intent(in) :: label integer :: it it = return_label_index(self, label) if (it /= 0) then call self%remove_entry_index(it) - else - return end if -end subroutine +end subroutine remove_entry_label -function get_n_entries(self) result(n) - class(double_dictionary_type) :: self +pure function get_n_entries(self) result(n) + class(double_dictionary_type), intent(in) :: self integer :: n n = self%n -end function +end function get_n_entries -subroutine copy(to, from) +pure subroutine copy(to, from) class(double_dictionary_type), intent(inout) :: to type(double_dictionary_type), intent(in) :: from integer :: n_entries, i to%n = from%n if (allocated(to%record)) deallocate(to%record) - if (from%get_n_entries() > 0) then + if (from%n > 0) then allocate(to%record(size(from%record))) - n_entries = from%get_n_entries() + n_entries = from%n do i = 1, n_entries to%record(i) = from%record(i) @@ -375,7 +321,26 @@ subroutine copy(to, from) end subroutine -subroutine copy_record(to, from) +pure subroutine move_record(from, to) + class(double_record), intent(inout) :: to + type(double_record), intent(inout) :: from + + call move_alloc(from%label, to%label) + if (allocated(to%array1)) deallocate(to%array1) + if (allocated(from%array1)) then + call move_alloc(from%array1, to%array1) + end if + if (allocated(to%array2)) deallocate(to%array2) + if (allocated(from%array2)) then + call move_alloc(from%array2, to%array2) + end if + if (allocated(to%array3)) deallocate(to%array3) + if (allocated(from%array3)) then + call move_alloc(from%array3, to%array3) + end if +end subroutine move_record + +pure subroutine copy_record(to, from) class(double_record), intent(inout) :: to type(double_record), intent(in) :: from @@ -386,7 +351,7 @@ subroutine copy_record(to, from) if (allocated(from%array3)) to%array3 = from%array3 end subroutine -subroutine concatenate_overwrite(self, dict2) +pure subroutine concatenate_overwrite(self, dict2) class(double_dictionary_type), intent(inout) :: self type(double_dictionary_type), intent(in) :: dict2 type(double_dictionary_type) :: tmp_dict @@ -395,14 +360,15 @@ subroutine concatenate_overwrite(self, dict2) end subroutine -function combine_dict(self, dict2) result(new_dict) +pure function combine_dict(self, dict2) result(new_dict) class(double_dictionary_type), intent(in) :: self type(double_dictionary_type), intent(in) :: dict2 type(double_dictionary_type) :: new_dict + integer :: it, i, n_entries new_dict = self associate(dict => dict2) - n_entries = dict%get_n_entries() + n_entries = dict%n do i = 1, n_entries call new_dict%push(dict%record(i)%label, it) new_dict%record(it) = dict%record(i) @@ -411,10 +377,11 @@ function combine_dict(self, dict2) result(new_dict) end function -subroutine update_1d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:) +pure subroutine update_1d(self, label, array) + class(double_dictionary_type), intent(inout) :: self + character(len=*), intent(in) :: label + real(wp), intent(in) :: array(:) + integer :: it it = return_label_index(self, label) @@ -430,10 +397,10 @@ subroutine update_1d(self, label, array) end if end subroutine -subroutine update_2d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:, :) +pure subroutine update_2d(self, label, array) + class(double_dictionary_type), intent(inout) :: self + character(len=*), intent(in) :: label + real(wp), intent(in) :: array(:, :) integer :: it it = return_label_index(self, label) @@ -449,10 +416,10 @@ subroutine update_2d(self, label, array) end if end subroutine -subroutine update_3d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:, :, :) +pure subroutine update_3d(self, label, array) + class(double_dictionary_type), intent(inout) :: self + character(len=*), intent(in) :: label + real(wp), intent(in) :: array(:, :, :) integer :: it it = return_label_index(self, label) @@ -463,116 +430,130 @@ subroutine update_3d(self, label, array) if (allocated(record%array3)) deallocate(record%array3) record%array3 = array end associate - else - return end if -end subroutine +end subroutine update_3d -subroutine get_label(self, index, label) - class(double_dictionary_type) :: self - integer :: index - character(len=:), allocatable :: label +!> Get the label of a given index +pure subroutine get_label(self, index, label) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Index of the label in the dictionary + integer, intent(in) :: index + !> Label to be returned + character(len=:), allocatable, intent(out) :: label if (index > self%n) return - if (allocated(label)) deallocate(label) label = self%record(index)%label -end subroutine +end subroutine get_label -subroutine get_1d_label(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp), allocatable :: array(:) +!> Get a 1D array from the dictionary +pure subroutine get_1d_label(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 1D array to be returned + real(wp), allocatable, intent(out) :: array(:) integer :: it it = return_label_index(self, label) if (it == 0) return call self%get_entry(it, array) -end subroutine +end subroutine get_1d_label -subroutine get_2d_label(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp), allocatable :: array(:,:) +!> Get a 2D array from the dictionary +pure subroutine get_2d_label(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 2D array to be returned + real(wp), allocatable, intent(out) :: array(:,:) integer :: it it = return_label_index(self, label) if (it == 0) return call self%get_entry(it, array) +end subroutine get_2d_label -end subroutine - -subroutine get_3d_label(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp), allocatable :: array(:, :, :) +!> Get a 3D array from the dictionary +pure subroutine get_3d_label(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 3D array to be returned + real(wp), allocatable, intent(out) :: array(:, :, :) integer :: it it = return_label_index(self, label) if (it == 0) return call self%get_entry(it, array) +end subroutine get_3d_label -end subroutine - -subroutine get_1d_index(self, index, array) - class(double_dictionary_type) :: self - integer :: index - real(wp), allocatable :: array(:) +!> Get a 1D array from the dictionary +pure subroutine get_1d_index(self, index, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Index of the array in the dictionary + integer, intent(in) :: index + !> 1D array to be returned + real(wp), allocatable, intent(out) :: array(:) if (index > self%n) return if (index <= 0) return associate(rec => self%record(index)) if (allocated(rec%array1)) then - if (allocated(array)) deallocate(array) - allocate(array(size(rec%array1, dim = 1))) array = rec%array1 - else - return end if end associate -end subroutine +end subroutine get_1d_index -subroutine get_2d_index(self, index, array) - class(double_dictionary_type) :: self - integer :: index - real(wp), allocatable :: array(:,:) +!> Get a 2D array from the dictionary +pure subroutine get_2d_index(self, index, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Index of the array in the dictionary + integer, intent(in) :: index + !> 2D array to be returned + real(wp), allocatable, intent(out) :: array(:,:) if (index > self%n) return if (index <= 0) return associate(rec => self%record(index)) if (allocated(rec%array2)) then - if (allocated(array)) deallocate(array) - allocate(array(size(rec%array2, dim = 1), size(rec%array2, dim = 2))) array = rec%array2 - else - return end if end associate -end subroutine +end subroutine get_2d_index -subroutine get_3d_index(self, index, array) - class(double_dictionary_type) :: self - integer :: index - real(wp), allocatable :: array(:, :, :) +!> Get a 3D array from the dictionary +pure subroutine get_3d_index(self, index, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(in) :: self + !> Index of the array in the dictionary + integer, intent(in) :: index + !> 3D array to be returned + real(wp), allocatable, intent(out) :: array(:, :, :) if (index > self%n) return if (index <= 0) return associate(rec => self%record(index)) if (allocated(rec%array3)) then - if (allocated(array)) deallocate(array) - allocate(array(size(rec%array3, dim = 1), size(rec%array3, dim = 2), size(rec%array3, dim=3))) array = rec%array3 - else - return end if end associate -end subroutine +end subroutine get_3d_index -subroutine push(self, label, it) +!> Add a new entry to the dictionary +pure subroutine push(self, label, it) + !> Instance of the double dictionary class(double_dictionary_type), intent(inout) :: self + !> Label for the array character(len=*), intent(in) :: label - + !> Index of the label in the dictionary integer, intent(out) :: it if (.not.allocated(self%record)) call resize(self%record) @@ -585,115 +566,148 @@ subroutine push(self, label, it) self%n = self%n + 1 it = self%n - self%record(it) = double_record(label=label) + self%record(it)%label = label end if end subroutine push -subroutine ini_label(self,label) - class(double_dictionary_type) :: self - character(len=*) :: label +!> Initialize a label in the dictionary +pure subroutine ini_label(self, label) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label integer :: it call self%push(label, it) -end subroutine +end subroutine ini_label -subroutine ini_1d(self, label, ndim1) - class(double_dictionary_type) :: self - character(len=*) :: label - integer :: ndim1 +!> Initialize a 1D array in the dictionary +pure subroutine ini_1d(self, label, ndim1) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> Dimensions of the 1D array + integer, intent(in) :: ndim1 integer :: it call self%push(label, it) associate(record => self%record(it)) - allocate(record%array1(ndim1), source = 0.0_wp) + allocate(record%array1(ndim1), source=0.0_wp) end associate -end subroutine - -subroutine ini_2d(self, label, ndim1, ndim2) +end subroutine ini_1d - class(double_dictionary_type) :: self - character(len=*) :: label - integer :: ndim1, ndim2 +!> Initialize a 2D array in the dictionary +pure subroutine ini_2d(self, label, ndim1, ndim2) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> Dimensions of the 2D array + integer, intent(in) :: ndim1, ndim2 integer :: it call self%push(label, it) associate(record => self%record(it)) - allocate(record%array2(ndim1, ndim2), source = 0.0_wp) + allocate(record%array2(ndim1, ndim2), source=0.0_wp) end associate -end subroutine +end subroutine ini_2d -subroutine ini_3d(self, label, ndim1, ndim2, ndim3) - class(double_dictionary_type) :: self - character(len=*) :: label - integer :: ndim1, ndim2, ndim3 +!> Initialize a 3D array in the dictionary +pure subroutine ini_3d(self, label, ndim1, ndim2, ndim3) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> Dimensions of the 3D array + integer, intent(in) :: ndim1, ndim2, ndim3 integer :: it call self%push(label, it) associate(record => self%record(it)) - allocate(record%array3(ndim1, ndim2, ndim3), source = 0.0_wp) + allocate(record%array3(ndim1, ndim2, ndim3), source=0.0_wp) end associate -end subroutine +end subroutine ini_3d -subroutine add_1d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:) +!> Add a 1D array to the dictionary +pure subroutine add_1d(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 1D array to be added + real(wp), intent(in) :: array(:) integer :: it call self%push(label, it) associate(record => self%record(it)) record%array1 = array + if (allocated(record%array2)) deallocate(record%array2) + if (allocated(record%array3)) deallocate(record%array3) end associate -end subroutine +end subroutine add_1d -subroutine add_2d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:,:) +!> Add a 2D array to the dictionary +pure subroutine add_2d(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 2D array to be added + real(wp), intent(in) :: array(:,:) integer :: it call self%push(label, it) associate(record => self%record(it)) record%array2 = array + if (allocated(record%array1)) deallocate(record%array1) + if (allocated(record%array3)) deallocate(record%array3) end associate -end subroutine +end subroutine add_2d -subroutine add_3d(self, label, array) - class(double_dictionary_type) :: self - character(len=*) :: label - real(wp) :: array(:, :, :) +!> Add a 3D array to the dictionary +pure subroutine add_3d(self, label, array) + !> Instance of the double dictionary + class(double_dictionary_type), intent(inout) :: self + !> Label for the array + character(len=*), intent(in) :: label + !> 3D array to be added + real(wp), intent(in) :: array(:, :, :) integer :: it call self%push(label, it) associate(record => self%record(it)) record%array3 = array + if (allocated(record%array1)) deallocate(record%array1) + if (allocated(record%array2)) deallocate(record%array2) end associate -end subroutine - +end subroutine add_3d -function return_label_index(self, label) result(it) +!> Find the index of a label in the dictionary +pure function return_label_index(self, label) result(it) class(double_dictionary_type), intent(in) :: self character(len=*), intent(in) :: label integer :: it it = 0 if (self%n <= 0) return it = find(self%record(:self%n), label) - if (it == 0) return - end function return_label_index - +!> Find the index of a label in the dictionary pure function find(record, label) result(pos) + !> List of double arrays type(double_record), intent(in) :: record(:) + !> Label to search for character(len=*), intent(in) :: label + !> Position of the label in the list integer :: pos integer :: i @@ -711,7 +725,7 @@ pure function find(record, label) result(pos) end function find - !> Reallocate list of double arrays +!> Reallocate list of double arrays pure subroutine resize(var, n) !> Instance of the array to be resized type(double_record), allocatable, intent(inout) :: var(:) @@ -719,7 +733,7 @@ pure subroutine resize(var, n) integer, intent(in), optional :: n type(double_record), allocatable :: tmp(:) - integer :: this_size, new_size + integer :: this_size, new_size, it integer, parameter :: initial_size = 20 if (allocated(var)) then @@ -739,10 +753,12 @@ pure subroutine resize(var, n) if (allocated(tmp)) then this_size = min(size(tmp, 1), size(var, 1)) - var(:this_size) = tmp(:this_size) + do it = 1, this_size + call move_record(tmp(it), var(it)) + end do deallocate(tmp) end if end subroutine resize -end module tblite_double_dictionary +end module tblite_double_dictionary \ No newline at end of file diff --git a/src/tblite/io/numpy/CMakeLists.txt b/src/tblite/io/numpy/CMakeLists.txt index 3e370a10..2fd1bfdf 100644 --- a/src/tblite/io/numpy/CMakeLists.txt +++ b/src/tblite/io/numpy/CMakeLists.txt @@ -25,6 +25,7 @@ list( "${dir}/save.f90" "${dir}/savez.f90" "${dir}/utils.f90" + "${dir}/zip.f90" ) set(srcs "${srcs}" PARENT_SCOPE) \ No newline at end of file diff --git a/src/tblite/io/numpy/crc32.f90 b/src/tblite/io/numpy/crc32.f90 index 276812ec..d90d6c1c 100644 --- a/src/tblite/io/numpy/crc32.f90 +++ b/src/tblite/io/numpy/crc32.f90 @@ -14,8 +14,8 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with tblite. If not, see . -!> @file tblite/io/numpy/load.f90 -!> Provides npy input routines +!> @file tblite/io/numpy/crc32.f90 +!> Provide crc32 hashing function !> Implementation of cyclic redundancy check hashing function, !> used to check the integrity of data in zip files. diff --git a/src/tblite/io/numpy/meson.build b/src/tblite/io/numpy/meson.build index 0d73d573..2c49e8e7 100644 --- a/src/tblite/io/numpy/meson.build +++ b/src/tblite/io/numpy/meson.build @@ -22,4 +22,5 @@ srcs += files( 'save.f90', 'savez.f90', 'utils.f90', + 'zip.f90', ) \ No newline at end of file diff --git a/src/tblite/io/numpy/savez.f90 b/src/tblite/io/numpy/savez.f90 index f79515b8..76d219a8 100644 --- a/src/tblite/io/numpy/savez.f90 +++ b/src/tblite/io/numpy/savez.f90 @@ -25,6 +25,7 @@ module tblite_io_numpy_savez & zip_global_sig, zip_local_sig, zip_footer_sig, zip_min_version use tblite_io_numpy_crc32, only : crc32_hash use tblite_io_numpy_save, only : npy_header + use tblite_io_numpy_zip, only : list_zip_file, zip_file use tblite_output_format, only : format_string implicit none private @@ -50,24 +51,23 @@ subroutine save_npz_i4_r1(filename, varname, array, iostat, iomsg) character(len=*), parameter :: vtype = type_i4 integer, parameter :: vsize = 4 - integer(i2) :: nrecs - integer(i4) :: global_header_offset, checksum + integer(i4) :: checksum, nbytes character(len=:), allocatable :: file_header, global_header, local_header, footer, path - integer(i4) :: nbytes logical :: exist integer :: io, stat character(len=:), allocatable :: msg + type(zip_file) :: zip path = varname // ".npy" - nrecs = 0_i2 - global_header_offset = 0_i4 - global_header = "" + zip%nrecs = 0_i2 + zip%global_header_offset = 0_i4 + zip%global_header = "" inquire(file=filename, exist=exist) open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) if (stat == 0 .and. exist) then - call parse_zip_footer(io, filename, nrecs, global_header, global_header_offset, stat, msg) + call list_zip_file(io, filename, zip, stat, msg) end if if (stat == 0) then @@ -79,10 +79,10 @@ subroutine save_npz_i4_r1(filename, varname, array, iostat, iomsg) nbytes = len(file_header) + size(array) * vsize local_header = get_local_header(path, checksum, nbytes) - global_header = global_header // get_global_header(path, local_header, global_header_offset) - footer = get_footer(nrecs + 1_i2, len(global_header), global_header_offset + len(local_header) + nbytes) + global_header = zip%global_header // get_global_header(path, local_header, zip%global_header_offset) + footer = get_footer(zip%nrecs + 1_i2, len(global_header), zip%global_header_offset + len(local_header) + nbytes) - write(io, pos=global_header_offset+1, iostat=stat) & + write(io, pos=zip%global_header_offset+1, iostat=stat) & & local_header, & & file_header, & & array, & @@ -105,24 +105,23 @@ subroutine save_npz_rdp_r1(filename, varname, array, iostat, iomsg) character(len=*), parameter :: vtype = type_rdp integer, parameter :: vsize = 8 - integer(i2) :: nrecs - integer(i4) :: global_header_offset, checksum + integer(i4) :: checksum, nbytes character(len=:), allocatable :: file_header, global_header, local_header, footer, path - integer(i4) :: nbytes logical :: exist integer :: io, stat character(len=:), allocatable :: msg + type(zip_file) :: zip path = varname // ".npy" - nrecs = 0_i2 - global_header_offset = 0_i4 - global_header = "" + zip%nrecs = 0_i2 + zip%global_header_offset = 0_i4 + zip%global_header = "" inquire(file=filename, exist=exist) open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) if (stat == 0 .and. exist) then - call parse_zip_footer(io, filename, nrecs, global_header, global_header_offset, stat, msg) + call list_zip_file(io, filename, zip, stat, msg) end if if (stat == 0) then @@ -134,10 +133,10 @@ subroutine save_npz_rdp_r1(filename, varname, array, iostat, iomsg) nbytes = len(file_header) + size(array) * vsize local_header = get_local_header(path, checksum, nbytes) - global_header = global_header // get_global_header(path, local_header, global_header_offset) - footer = get_footer(nrecs + 1_i2, len(global_header), global_header_offset + len(local_header) + nbytes) + global_header = zip%global_header // get_global_header(path, local_header, zip%global_header_offset) + footer = get_footer(zip%nrecs + 1_i2, len(global_header), zip%global_header_offset + len(local_header) + nbytes) - write(io, pos=global_header_offset+1, iostat=stat) & + write(io, pos=zip%global_header_offset+1, iostat=stat) & & local_header, & & file_header, & & array, & @@ -160,24 +159,23 @@ subroutine save_npz_rdp_r2(filename, varname, array, iostat, iomsg) character(len=*), parameter :: vtype = type_rdp integer, parameter :: vsize = 8 - integer(i2) :: nrecs - integer(i4) :: global_header_offset, checksum + integer(i4) :: checksum, nbytes character(len=:), allocatable :: file_header, global_header, local_header, footer, path - integer(i4) :: nbytes logical :: exist integer :: io, stat character(len=:), allocatable :: msg + type(zip_file) :: zip path = varname // ".npy" - nrecs = 0_i2 - global_header_offset = 0_i4 - global_header = "" + zip%nrecs = 0_i2 + zip%global_header_offset = 0_i4 + zip%global_header = "" inquire(file=filename, exist=exist) open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) if (stat == 0 .and. exist) then - call parse_zip_footer(io, filename, nrecs, global_header, global_header_offset, stat, msg) + call list_zip_file(io, filename, zip, stat, msg) end if if (stat == 0) then @@ -189,10 +187,10 @@ subroutine save_npz_rdp_r2(filename, varname, array, iostat, iomsg) nbytes = len(file_header) + size(array) * vsize local_header = get_local_header(path, checksum, nbytes) - global_header = global_header // get_global_header(path, local_header, global_header_offset) - footer = get_footer(nrecs + 1_i2, len(global_header), global_header_offset + len(local_header) + nbytes) + global_header = zip%global_header // get_global_header(path, local_header, zip%global_header_offset) + footer = get_footer(zip%nrecs + 1_i2, len(global_header), zip%global_header_offset + len(local_header) + nbytes) - write(io, pos=global_header_offset+1, iostat=stat) & + write(io, pos=zip%global_header_offset+1, iostat=stat) & & local_header, & & file_header, & & array, & @@ -215,24 +213,23 @@ subroutine save_npz_rdp_r3(filename, varname, array, iostat, iomsg) character(len=*), parameter :: vtype = type_rdp integer, parameter :: vsize = 8 - integer(i2) :: nrecs - integer(i4) :: global_header_offset, checksum + integer(i4) :: checksum, nbytes character(len=:), allocatable :: file_header, global_header, local_header, footer, path - integer(i4) :: nbytes logical :: exist integer :: io, stat character(len=:), allocatable :: msg + type(zip_file) :: zip path = varname // ".npy" - nrecs = 0_i2 - global_header_offset = 0_i4 - global_header = "" + zip%nrecs = 0_i2 + zip%global_header_offset = 0_i4 + zip%global_header = "" inquire(file=filename, exist=exist) open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) if (stat == 0 .and. exist) then - call parse_zip_footer(io, filename, nrecs, global_header, global_header_offset, stat, msg) + call list_zip_file(io, filename, zip, stat, msg) end if if (stat == 0) then @@ -244,10 +241,10 @@ subroutine save_npz_rdp_r3(filename, varname, array, iostat, iomsg) nbytes = len(file_header) + size(array) * vsize local_header = get_local_header(path, checksum, nbytes) - global_header = global_header // get_global_header(path, local_header, global_header_offset) - footer = get_footer(nrecs + 1_i2, len(global_header), global_header_offset + len(local_header) + nbytes) + global_header = zip%global_header // get_global_header(path, local_header, zip%global_header_offset) + footer = get_footer(zip%nrecs + 1_i2, len(global_header), zip%global_header_offset + len(local_header) + nbytes) - write(io, pos=global_header_offset+1, iostat=stat) & + write(io, pos=zip%global_header_offset+1, iostat=stat) & & local_header, & & file_header, & & array, & @@ -339,146 +336,6 @@ pure function get_footer(nrecs, global_header_size, global_header_offset) result end function get_footer -subroutine parse_zip_footer(io, filename, nrecs, global_header, global_header_offset, stat, msg) - integer, intent(in) :: io - character(len=*), intent(in) :: filename - integer(i2), intent(out) :: nrecs - character(len=:), allocatable, intent(out) :: global_header - integer(i4), intent(out) :: global_header_offset - integer, intent(out) :: stat - character(len=:), allocatable :: msg - - integer(i2) :: path_size, extra_field_size, comment_size - integer(i2) :: disk_no, disk_start, nrecs_on_disk - integer(i4) :: nbytes_compressed, global_header_size - character(len=512) :: errmsg - integer :: res, length, pos - integer(i4) :: header_sig - character(len=:), allocatable :: path - - stat = 0 - pos = 1 - read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig - do while(stat == 0 .and. is_local_header(header_sig)) - - if (stat == 0) & - read(io, pos=pos+18, iostat=stat, iomsg=errmsg) nbytes_compressed - if (stat == 0) & - read(io, pos=pos+26, iostat=stat, iomsg=errmsg) path_size - if (stat == 0) & - read(io, pos=pos+28, iostat=stat, iomsg=errmsg) extra_field_size - - if (stat == 0) then - if (allocated(path)) deallocate(path) - allocate(character(len=path_size) :: path, stat=stat) - end if - if (stat == 0) & - read(io, pos=pos+30, iostat=stat, iomsg=errmsg) path - - pos = pos + 30 + path_size + extra_field_size + nbytes_compressed - read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig - end do - if (stat /= 0) then - msg = "Failed to read local header block for '"//filename//"'" - if (len_trim(errmsg) > 0) & - msg = msg // " ("//trim(errmsg)//")" - return - end if - - if (.not.is_global_header(header_sig)) then - stat = 400 - msg = "Expected global header signature for '"//filename//"' got "// & - & format_string(header_sig, '(z0.8)') - return - end if - - ! global_header_offset = pos - 1 - do while(stat == 0 .and. is_global_header(header_sig)) - - if (stat == 0) & - read(io, pos=pos+28, iostat=stat, iomsg=errmsg) path_size - if (stat == 0) & - read(io, pos=pos+30, iostat=stat, iomsg=errmsg) extra_field_size - if (stat == 0) & - read(io, pos=pos+32, iostat=stat, iomsg=errmsg) comment_size - - if (stat == 0) then - if (allocated(path)) deallocate(path) - allocate(character(len=path_size) :: path, stat=stat) - end if - if (stat == 0) & - read(io, pos=pos+46, iostat=stat, iomsg=errmsg) path - - pos = pos + 46 + path_size + extra_field_size + comment_size - read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig - end do - if (stat /= 0) then - msg = "Failed to read global header block for '"//filename//"'" - if (len_trim(errmsg) > 0) & - msg = msg // " ("//trim(errmsg)//")" - return - end if - if (.not.is_footer_header(header_sig)) then - stat = 401 - msg = "Expected footer signature for '"//filename//"' got "// & - & format_string(header_sig, '(z0.8)') - return - end if - ! global_header_size = pos - global_header_offset + 1 - - read(io, pos=pos+4, iostat=stat, iomsg=errmsg) & - & disk_no, disk_start, nrecs_on_disk, nrecs, & - & global_header_size, global_header_offset, comment_size - - if (stat == 0) & - allocate(character(len=global_header_size) :: global_header, stat=stat) - if (stat == 0) & - read(io, iostat=stat, pos=global_header_offset+1) global_header - - if (stat /= 0) then - msg = "Failed to read footer for '"//filename//"'" - if (len_trim(errmsg) > 0) & - msg = msg // " ("//trim(errmsg)//")" - return - end if - - if (disk_no /= 0) then - stat = 402 - msg = "Cannot read zip file with disk_no != 0" - end if - - if (disk_start /= 0) then - stat = 403 - msg = "Cannot read zip file with disk_start != 0" - end if - - if (nrecs_on_disk /= nrecs) then - stat = 404 - msg = "Cannot read zip file with nrecs_on_disk != nrecs" - end if -end subroutine parse_zip_footer - -pure function is_local_header(header_sig) result(is_local) - integer(i4), intent(in) :: header_sig - logical :: is_local - - is_local = header_sig == zip_local_sig -end function is_local_header - -pure function is_global_header(header_sig) result(is_global) - integer(i4), intent(in) :: header_sig - logical :: is_global - - is_global = header_sig == zip_global_sig -end function is_global_header - -pure function is_footer_header(header_sig) result(is_footer) - integer(i4), intent(in) :: header_sig - logical :: is_footer - - is_footer = header_sig == zip_footer_sig -end function is_footer_header - !> Handle iostat and iomsg of write operation subroutine handle_iostat(stat, filename, iostat, iomsg) !> Error status of loading, zero on success diff --git a/src/tblite/io/numpy/zip.f90 b/src/tblite/io/numpy/zip.f90 new file mode 100644 index 00000000..3573a947 --- /dev/null +++ b/src/tblite/io/numpy/zip.f90 @@ -0,0 +1,191 @@ +! This file is part of tblite. +! SPDX-Identifier: LGPL-3.0-or-later +! +! tblite is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! tblite is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with tblite. If not, see . + +!> @file tblite/io/numpy/zip.f90 +!> Provides zip implementation for numpy input/output. + +!> Provide routines for handling zip files +module tblite_io_numpy_zip + use mctc_env, only : i2, i4 + use tblite_io_numpy_constants, only : & + & zip_global_sig, zip_local_sig, zip_footer_sig, zip_min_version + use tblite_output_format, only : format_string + implicit none + private + + public :: list_zip_file, zip_file + + type :: zip_record + character(len=:), allocatable :: path + end type zip_record + + type :: zip_file + type(zip_record), allocatable :: records(:) + character(len=:), allocatable :: global_header + integer(i4) :: global_header_offset = 0_i4 + integer(i2) :: nrecs = 0_i2 + end type zip_file + + +contains + +subroutine list_zip_file(io, filename, zip, stat, msg) + integer, intent(in) :: io + character(len=*), intent(in) :: filename + type(zip_file), intent(out) :: zip + integer, intent(out) :: stat + character(len=:), allocatable :: msg + + integer :: irec + integer(i2) :: path_size, extra_field_size, comment_size + integer(i2) :: disk_no, disk_start, nrecs_on_disk + integer(i4) :: nbytes_compressed, global_header_size + character(len=512) :: errmsg + integer :: res, length, pos + integer(i4) :: header_sig + character(len=:), allocatable :: path + + stat = 0 + irec = 0 + pos = 1 + read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig + do while(stat == 0 .and. is_local_header(header_sig)) + irec = irec + 1 + + if (stat == 0) & + read(io, pos=pos+18, iostat=stat, iomsg=errmsg) nbytes_compressed + if (stat == 0) & + read(io, pos=pos+26, iostat=stat, iomsg=errmsg) path_size + if (stat == 0) & + read(io, pos=pos+28, iostat=stat, iomsg=errmsg) extra_field_size + + if (stat == 0) then + if (allocated(path)) deallocate(path) + allocate(character(len=path_size) :: path, stat=stat) + end if + if (stat == 0) & + read(io, pos=pos+30, iostat=stat, iomsg=errmsg) path + + pos = pos + 30 + path_size + extra_field_size + nbytes_compressed + read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig + end do + if (stat /= 0) then + msg = "Failed to read local header block for '"//filename//"'" + if (len_trim(errmsg) > 0) & + msg = msg // " ("//trim(errmsg)//")" + return + end if + + if (.not.is_global_header(header_sig)) then + stat = 400 + msg = "Expected global header signature for '"//filename//"' got "// & + & format_string(header_sig, '(z0.8)') + return + end if + + allocate(zip%records(irec)) + + irec = 0 + ! global_header_offset = pos - 1 + do while(stat == 0 .and. is_global_header(header_sig)) + irec = irec + 1 + + if (stat == 0) & + read(io, pos=pos+28, iostat=stat, iomsg=errmsg) path_size + if (stat == 0) & + read(io, pos=pos+30, iostat=stat, iomsg=errmsg) extra_field_size + if (stat == 0) & + read(io, pos=pos+32, iostat=stat, iomsg=errmsg) comment_size + + if (stat == 0) then + if (allocated(path)) deallocate(path) + allocate(character(len=path_size) :: path, stat=stat) + end if + if (stat == 0) & + read(io, pos=pos+46, iostat=stat, iomsg=errmsg) path + + zip%records(irec)%path = path + + pos = pos + 46 + path_size + extra_field_size + comment_size + read(io, pos=pos, iostat=stat, iomsg=errmsg) header_sig + end do + if (stat /= 0) then + msg = "Failed to read global header block for '"//filename//"'" + if (len_trim(errmsg) > 0) & + msg = msg // " ("//trim(errmsg)//")" + return + end if + if (.not.is_footer_header(header_sig)) then + stat = 401 + msg = "Expected footer signature for '"//filename//"' got "// & + & format_string(header_sig, '(z0.8)') + return + end if + ! global_header_size = pos - global_header_offset + 1 + + read(io, pos=pos+4, iostat=stat, iomsg=errmsg) & + & disk_no, disk_start, nrecs_on_disk, zip%nrecs, & + & global_header_size, zip%global_header_offset, comment_size + + if (stat == 0) & + allocate(character(len=global_header_size) :: zip%global_header, stat=stat) + if (stat == 0) & + read(io, iostat=stat, pos=zip%global_header_offset+1) zip%global_header + + if (stat /= 0) then + msg = "Failed to read footer for '"//filename//"'" + if (len_trim(errmsg) > 0) & + msg = msg // " ("//trim(errmsg)//")" + return + end if + + if (disk_no /= 0) then + stat = 402 + msg = "Cannot read zip file with disk_no != 0" + end if + + if (disk_start /= 0) then + stat = 403 + msg = "Cannot read zip file with disk_start != 0" + end if + + if (nrecs_on_disk /= zip%nrecs) then + stat = 404 + msg = "Cannot read zip file with nrecs_on_disk != nrecs" + end if +end subroutine list_zip_file + +pure function is_local_header(header_sig) result(is_local) + integer(i4), intent(in) :: header_sig + logical :: is_local + + is_local = header_sig == zip_local_sig +end function is_local_header + +pure function is_global_header(header_sig) result(is_global) + integer(i4), intent(in) :: header_sig + logical :: is_global + + is_global = header_sig == zip_global_sig +end function is_global_header + +pure function is_footer_header(header_sig) result(is_footer) + integer(i4), intent(in) :: header_sig + logical :: is_footer + + is_footer = header_sig == zip_footer_sig +end function is_footer_header +end module tblite_io_numpy_zip \ No newline at end of file diff --git a/src/tblite/xtb/singlepoint.f90 b/src/tblite/xtb/singlepoint.f90 index 11c5989a..c11bd83c 100644 --- a/src/tblite/xtb/singlepoint.f90 +++ b/src/tblite/xtb/singlepoint.f90 @@ -345,7 +345,6 @@ subroutine xtb_singlepoint(ctx, mol, calc, wfn, accuracy, energy, gradient, sigm allocate(results%dict) if (present(post_process)) then call post_process%pack_res(mol, results) - if (prlevel > 1) call results%dict%dump("post_processing.toml", error) end if end if diff --git a/test/unit/test_double_dictionary.f90 b/test/unit/test_double_dictionary.f90 index da389147..c7c48dbe 100644 --- a/test/unit/test_double_dictionary.f90 +++ b/test/unit/test_double_dictionary.f90 @@ -3,6 +3,7 @@ module test_double_dictionary use mctc_env_testing, only : new_unittest, unittest_type, error_type, check, & & test_failed use tblite_double_dictionary, only : double_dictionary_type + use tblite_io_numpy, only : save_npz implicit none private @@ -28,8 +29,8 @@ subroutine collect_double_dictionary(testsuite) new_unittest("compare index and label lookup", test_equivalence_index_label_lookup), & new_unittest("initialize labels", test_initialize_labels), & new_unittest("update entries", test_update_entries_label), & - new_unittest("read in toml-structured file", test_read_in_toml), & - new_unittest("write toml-structured output and read in again", test_write_read_toml), & + new_unittest("read in npz file", test_read_in_npz), & + new_unittest("write npz output and read in again", test_write_read_npz), & new_unittest("check equal operator", test_equal_operator), & new_unittest("check equal operator for two different dicts", test_equal_operator_different_dict) & ] @@ -497,6 +498,73 @@ subroutine test_addition_operator(error) end subroutine +subroutine test_read_in_npz(error) + type(error_type), allocatable, intent(out) :: error + type(double_dictionary_type) :: dict1, dict2 + character(len=*), parameter :: filename = ".read-ddict.npz" + + real(wp), allocatable :: input1(:), input2(:, :), input3(:, :, :) + real(wp), allocatable :: array1(:), array2(:, :), array3(:, :, :) + allocate(input1(4), source = 1.0_wp) + allocate(input2(4, 6), source = 2.0_wp) + allocate(input3(4, 6, 9), source = 0.0_wp) + + call delete_file(filename) + call save_npz(filename, "tblite0_r1_a", input1) + call save_npz(filename, "tblite0_r2_b", input2) + call save_npz(filename, "tblite0_r3_c", input3) + + call dict2%load(filename, error) + call delete_file(filename) + if (allocated(error)) return + + call dict2%get_entry("a", array1) + call check(error, allocated(array1), "First entry not allocated") + if (allocated(error)) return + call dict2%get_entry("b", array2) + call check(error, allocated(array2), "Second entry not allocated") + if (allocated(error)) return + call dict2%get_entry("c", array3) + call check(error, allocated(array3), "Third entry not allocated") + if (allocated(error)) return + + call check(error, sum(array1 - input1), 0.0_wp) + if (allocated(error)) return + call check(error, sum(array2 - input2), 0.0_wp) + if (allocated(error)) return + call check(error, sum(array3 - input3), 0.0_wp) + if (allocated(error)) return +end subroutine test_read_in_npz + + +subroutine test_write_read_npz(error) + type(error_type), allocatable, intent(out) :: error + type(double_dictionary_type) :: dict1, dict2 + character(len=*), parameter :: filename = ".read-write-ddict.npz" + real(wp), allocatable :: array1(:), array2(:, :), array3(:, :, :) + + call fill_test_dict(dict1) + + call dict1%dump(filename, error) + if (allocated(error)) return + + call dict2%load(filename, error) + call delete_file(filename) + if (allocated(error)) return + + call dict2%get_entry("test1", array1) + call check(error, allocated(array1), "First entry not allocated") + if (allocated(error)) return + call dict2%get_entry("test2", array2) + call check(error, allocated(array2), "Second entry not allocated") + if (allocated(error)) return + call dict2%get_entry("test3", array3) + call check(error, allocated(array3), "Third entry not allocated") + if (allocated(error)) return + + call check(error, dict1 == dict2) +end subroutine test_write_read_npz + subroutine test_update_entries_label(error) type(error_type), allocatable, intent(out) :: error @@ -554,79 +622,6 @@ subroutine test_update_entries_label(error) if (allocated(error)) return end subroutine -subroutine test_write_read_toml(error) - type(error_type), allocatable, intent(out) :: error - type(double_dictionary_type) :: dict, dict1 - integer :: io = 42 - - call fill_test_dict(dict) - open(newunit=io, status="scratch") - call dict%dump(io, error) - rewind io - dict = double_dictionary_type(record=null()) - call fill_test_dict_1d_array(dict) - call dict1%load(2, error) - if (.not.(allocated(error))) then - allocate(error) - error%message = "Non-existent unit was opened" - return - end if - deallocate(error) - call dict1%load(io, error) - close(io) - call check(error, (dict == dict1)) - - call dict%dump("test.toml", error) - dict1 = double_dictionary_type(record=null()) - call dict1%load("test.tom", error) - if (.not.(allocated(error))) then - allocate(error) - error%message = "Non-existent file was opened" - return - end if - deallocate(error) - call dict1%load("test.toml", error) - call delete_file("test.toml") - - call check(error, (dict == dict1)) - -end subroutine - -subroutine test_read_in_toml(error) - type(error_type), allocatable, intent(out) :: error - type(double_dictionary_type) :: dict - integer :: io - - open(newunit=io, status="scratch") - write(io, '(a)') & - "levels = [ -1.3970922000000002E+01, -1.0063292000000001E+01 ]", & - "slater = [ 2.0964320000000001E+00, 1.8000000000000000E+00 ]", & - "" - rewind io - - call dict%load(io, error) - close(io) - - call check(error, (dict%get_n_entries() == 2)) - - open(newunit=io, status="scratch") - write(io, '(a)') & - "levels = [ -1.3970922000000002E+01, -1.0063292000000001E+01 ]", & - "slater = [ 2.0964320000000001E+00, 1.8000000000000000E+00 ", & - "" - rewind io - - call dict%load(io, error) - close(io) - if (.not.(allocated(error))) then - allocate(error) - error%message = "Faulty toml unit was parsed" - return - end if - deallocate(error) - -end subroutine - subroutine test_equal_operator(error) type(error_type), allocatable, intent(out) :: error type(double_dictionary_type) :: dict, dict_ diff --git a/test/unit/test_xtbml.f90 b/test/unit/test_xtbml.f90 index 8762c44d..4fdaf78a 100644 --- a/test/unit/test_xtbml.f90 +++ b/test/unit/test_xtbml.f90 @@ -13,6 +13,7 @@ module test_xtbml use tblite_data_spin, only : get_spin_constant use tblite_double_dictionary, only : double_dictionary_type use tblite_integral_type, only : integral_type + use tblite_io_numpy, only : save_npz use tblite_param_post_processing, only : post_processing_param_list use tblite_param_serde, only : serde_record use tblite_param_xtbml_features, only : xtbml_features_record @@ -696,41 +697,35 @@ subroutine test_orbital_energy_ref(error) !> Error handling integer,parameter :: nat = 3 type(error_type), allocatable, intent(out) :: error - real(wp), parameter :: xyz(3, nat) = reshape((/0.00000,0.00000,1.0000000,& - &0.0000,0.00000,-1.0000000,& - &0.000000,0.000000,0.000000/),shape=(/3,nat/)) + real(wp), parameter :: xyz(3, nat) = reshape([& + & 0.000000_wp, 0.000000_wp, 1.0000000_wp,& + & 0.000000_wp, 0.000000_wp, -1.0000000_wp,& + & 0.000000_wp, 0.000000_wp, 0.0000000_wp], shape=shape(xyz)) integer, parameter :: num(nat) = (/8,8,6/) type(results_type) :: res, res_ real(wp) :: energy = 0.0_wp real(wp), allocatable :: xtbml(:,:),xtbml_rot(:,:),xtbml_trans(:,:) integer :: i - integer :: io - - open(newunit=io, status="scratch") - write(io, '(a)') & - "response_alpha = [ 0.1337835451437173, 0.1337835451437157, 1.6307951377273491 ]", & - "gap_alpha = [ 0.0095881103175756, 0.0095881103175756, 0.0011085237545836 ]", & - "chem_pot_alpha = [ -0.9622382191942670, -0.9622382191942677, -0.8758044318869692 ]", & - "HOAO_alpha = [ -0.9670322743530548, -0.9670322743530555, -0.8763586937642610 ]", & - "LUAO_alpha = [ -0.9574441640354792, -0.9574441640354798, -0.8752501700096774 ]", & - "response_beta = [ 1.2449917888491431, 1.2449917888491182, 0.0059865937960335 ]", & - "gap_beta = [ 0.0806953599842188, 0.0806953599842184, 19.4997052573882037 ]", & - "chem_pot_beta = [ -15.1513328627176183, -15.1513328627176254, -8.6913073970860637 ]", & - "HOAO_beta = [ -15.1916805427097277, -15.1916805427097348, -18.4411600257801638 ]", & - "LUAO_beta = [ -15.1109851827255088, -15.1109851827255159, 1.0585452316080382 ]", & - "ext_gap_alpha = [ 0.0054951699014186, 0.0054951699014186, 0.0067615787201293 ]", & - "ext_chem_pot_alpha = [ -0.9205182203755125, -0.9205182203755127, -0.9334269287974830 ]", & - "ext_HOAO_alpha = [ -0.9232658053262217, -0.9232658053262219, -0.9368077181575476 ]", & - "ext_LUAO_alpha = [ -0.9177706354248032, -0.9177706354248034, -0.9300461394374184 ]", & - "ext_gap_beta = [ 9.4593893126801003, 9.4593893126801003, 6.5604665198186094 ]", & - "ext_chem_pot_beta = [ -12.9537156502782160, -12.9537156502782196, -13.9314158798429979 ]", & - "ext_HOAO_beta = [ -17.6834103066182671, -17.6834103066182706, -17.2116491397523035 ]", & - "ext_LUAO_beta = [ -8.2240209939381650, -8.2240209939381685, -10.6511826199336923 ]", & - "" - rewind io - call dict_ref%load(io, error) - close(io) + + call dict_ref%add_entry("response_alpha", [ 0.1337835451437173_wp, 0.1337835451437157_wp, 1.6307951377273491_wp]) + call dict_ref%add_entry("gap_alpha", [ 0.0095881103175756_wp, 0.0095881103175756_wp, 0.0011085237545836_wp]) + call dict_ref%add_entry("chem_pot_alpha", [ -0.9622382191942670_wp, -0.9622382191942677_wp, -0.8758044318869692_wp]) + call dict_ref%add_entry("HOAO_alpha", [ -0.9670322743530548_wp, -0.9670322743530555_wp, -0.8763586937642610_wp]) + call dict_ref%add_entry("LUAO_alpha", [ -0.9574441640354792_wp, -0.9574441640354798_wp, -0.8752501700096774_wp]) + call dict_ref%add_entry("response_beta", [ 1.2449917888491431_wp, 1.2449917888491182_wp, 0.0059865937960335_wp]) + call dict_ref%add_entry("gap_beta", [ 0.0806953599842188_wp, 0.0806953599842184_wp, 19.4997052573882037_wp]) + call dict_ref%add_entry("chem_pot_beta", [ -15.1513328627176183_wp, -15.1513328627176254_wp, -8.6913073970860637_wp]) + call dict_ref%add_entry("HOAO_beta", [ -15.1916805427097277_wp, -15.1916805427097348_wp, -18.4411600257801638_wp]) + call dict_ref%add_entry("LUAO_beta", [ -15.1109851827255088_wp, -15.1109851827255159_wp, 1.0585452316080382_wp]) + call dict_ref%add_entry("ext_gap_alpha", [ 0.0054951699014186_wp, 0.0054951699014186_wp, 0.0067615787201293_wp]) + call dict_ref%add_entry("ext_chem_pot_alpha", [ -0.9205182203755125_wp, -0.9205182203755127_wp, -0.9334269287974830_wp]) + call dict_ref%add_entry("ext_HOAO_alpha", [ -0.9232658053262217_wp, -0.9232658053262219_wp, -0.9368077181575476_wp]) + call dict_ref%add_entry("ext_LUAO_alpha", [ -0.9177706354248032_wp, -0.9177706354248034_wp, -0.9300461394374184_wp]) + call dict_ref%add_entry("ext_gap_beta", [ 9.4593893126801003_wp, 9.4593893126801003_wp, 6.5604665198186094_wp]) + call dict_ref%add_entry("ext_chem_pot_beta", [ -12.9537156502782160_wp, -12.9537156502782196_wp, -13.9314158798429979_wp]) + call dict_ref%add_entry("ext_HOAO_beta", [ -17.6834103066182671_wp, -17.6834103066182706_wp, -17.2116491397523035_wp]) + call dict_ref%add_entry("ext_LUAO_beta", [ -8.2240209939381650_wp, -8.2240209939381685_wp, -10.6511826199336923_wp]) call new(mol,num,xyz*aatoau,uhf=2,charge=0.0_wp) call new_gfn2_calculator(calc,mol, error) @@ -1080,10 +1075,10 @@ subroutine test_energy_sum_up_gfn2(error) call res%dict%get_entry("E_tot", tmp_array) sum_energy = sum_energy + sum(tmp_array) - print'(3es21.14)', abs(sum_energy-energy) + ! print'(3es21.14)', abs(sum_energy-energy) if (abs(sum_energy-energy) > thr) then call test_failed(error, "GFN2: Energy features don't add up to total energy.") - print'(3es21.14)', abs(sum_energy-energy) + ! print'(3es21.14)', abs(sum_energy-energy) end if @@ -1097,10 +1092,10 @@ subroutine test_energy_sum_up_gfn2(error) sum_energy = sum_energy + sum(tmp_array) deallocate(tmp_array) end do - print'(3es21.14)', abs(sum_energy-energy) + ! print'(3es21.14)', abs(sum_energy-energy) if (abs(sum_energy-energy) > thr) then call test_failed(error, "GFN2: Energy features don't add up to total energy.") - print'(3es21.14)', abs(sum_energy-energy) + ! print'(3es21.14)', abs(sum_energy-energy) end if end subroutine test_energy_sum_up_gfn2