diff --git a/CMakeLists.txt b/CMakeLists.txt index f476b20107..ff79c2085b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -131,6 +131,7 @@ option(CP2K_USE_GREENX "Enable GreenX support" ${CP2K_USE_EVERYTHING}) option(CP2K_USE_HDF5 "Enable HDF5 support" ${CP2K_USE_EVERYTHING}) option(CP2K_USE_LIBINT2 "Enable Libint2 support" ${CP2K_USE_EVERYTHING}) option(CP2K_USE_LIBTORCH "Enable LibTorch support" ${CP2K_USE_EVERYTHING}) +option(CP2K_USE_OPENPMD "Enable support for I/O via openPMD" OFF) option(CP2K_USE_LIBXC "Enable LibXC support" ${CP2K_USE_EVERYTHING}) option(CP2K_USE_MPI "Enable MPI support" ${CP2K_USE_EVERYTHING}) option(CP2K_USE_PEXSI "Enable PEXSI support" ${CP2K_USE_EVERYTHING}) @@ -821,6 +822,16 @@ if(CP2K_USE_MIMIC) find_package(MiMiC REQUIRED) endif() +if(CP2K_USE_OPENPMD) + find_package(openPMD 0.16.1 REQUIRED) + if(NOT openPMD_HAVE_MPI) + message( + FATAL_ERROR + "Found serial build of openPMD at ${openPMD_DIR}. Provide either an MPI build or deactivate openPMD." + ) + endif() +endif() + if(CP2K_USE_MPI_F08 AND NOT MPI_Fortran_HAVE_F08_MODULE) message( FATAL_ERROR @@ -983,6 +994,10 @@ if(CP2K_USE_HDF5) " - libraries: ${HDF5_LIBRARIES}\n\n") endif() +if(CP2K_USE_OPENPMD) + message(" - openPMD\n" " - libraries : ${openPMD_DIR}\n\n") +endif() + if(CP2K_USE_FFTW3) message(" - FFTW3\n" " - include directories: ${CP2K_FFTW3_INCLUDE_DIRS}\n" @@ -1158,6 +1173,10 @@ if(NOT CP2K_USE_HDF5) message(" - HDF5") endif() +if(NOT CP2K_USE_OPENPMD) + message(" - openPMD") +endif() + if(${CP2K_USE_ACCEL} MATCHES "NONE") message(" - GPU acceleration is disabled") endif() diff --git a/cmake/cmake_cp2k.sh b/cmake/cmake_cp2k.sh index 1edf9a7598..a059024ea0 100644 --- a/cmake/cmake_cp2k.sh +++ b/cmake/cmake_cp2k.sh @@ -51,6 +51,7 @@ if [[ "${PROFILE}" == "spack" ]] && [[ "${VERSION}" == "psmp" ]]; then -DCP2K_USE_DLAF=OFF \ -DCP2K_USE_LIBXSMM=OFF \ -DCP2K_USE_TBLITE=OFF \ + -DCP2K_USE_OPENPMD=ON \ -Werror=dev \ .. |& tee ./cmake.log CMAKE_EXIT_CODE=$? diff --git a/docs/technologies/libraries.md b/docs/technologies/libraries.md index 4730a92eed..d3706acfe8 100644 --- a/docs/technologies/libraries.md +++ b/docs/technologies/libraries.md @@ -229,6 +229,18 @@ greenX - Open-source file format and library. Support for greenX can be enabled - For more information see - Pass `-DCP2K_USE_TBLITE=ON` to CMake. +## openPMD (structured output) + +openPMD - Open-source data standard and library. Support for openPMD can be enabled in CMake via +`-DCP2K_USE_OPENPMD=ON`. CMake is the only supported way of enabling openPMD, use of `-D__OPENPMD` +as part of DFLAGS may or may not work. + +- openPMD-api may be downloaded from , a equal to or + greater than 0.16.1 is required. +- For more information see . +- The version of openPMD-api, determined by OPENPMDAPI_VERSION_GE, must be 0.16.1 or greater. +- openPMD-api must be built against MPI, determined by openPMD_HAVE_MPI. + ## HDF5 - Pass `-DCP2K_USE_HDF5=ON` to CMake to enable HDF5 support. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 154ffe9b0d..15152a09ca 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1561,6 +1561,10 @@ endif() list(APPEND CP2K_SRCS_C ${CP2K_DBM_SRCS_C} ${CP2K_GRID_SRCS_C}) +list(APPEND CP2K_SRCS_F openPMD.F input/cp_output_handling_openpmd.F + cp_realspace_grid_openpmd.F pw/realspace_grid_openpmd.F) +list(APPEND CP2K_SRCS_CPP openPMD.cpp) + # ############################################################################## # check that the files registered in CMakeLists.txt are all present and return # an error otherwise. cmake only runs this test nothing else below @@ -1659,15 +1663,12 @@ foreach(cp2k_prog ${CP2K_PROGS_F}) ${cp2k_prog} PROPERTIES COMPILE_DEFINITIONS __SHORT_FILE__="${cp2k_prog}") endforeach() -if(CP2K_USE_LIBTORCH) - add_library(cp2k +if(BUILD_SHARED_LIBS) + add_library(cp2k SHARED "${CP2K_SRCS};${CP2K_SRCS_C};${CP2K_SRCS_GPU};${CP2K_SRCS_CPP}") else() - if(BUILD_SHARED_LIBS) - add_library(cp2k SHARED "${CP2K_SRCS};${CP2K_SRCS_C};${CP2K_SRCS_GPU}") - else() - add_library(cp2k STATIC "${CP2K_SRCS};${CP2K_SRCS_C};${CP2K_SRCS_GPU}") - endif() + add_library(cp2k STATIC + "${CP2K_SRCS};${CP2K_SRCS_C};${CP2K_SRCS_GPU};${CP2K_SRCS_CPP}") endif() target_compile_features(cp2k PUBLIC cuda_std_11) @@ -1751,6 +1752,7 @@ target_link_libraries( cp2k_linalg_libs # Torch includes some LAPACK routines, but with floating-point exceptions. $<$:${TORCH_LIBRARIES}> + $<$:openPMD::openPMD> $<$,$>:gcov>) string(TIMESTAMP CP2K_TIMESTAMP "%Y-%m-%d %H:%M:%S") @@ -1790,6 +1792,7 @@ target_compile_definitions( $<$:__DLAF> $<$:__LIBXC> $<$:__HDF5> + $<$:__OPENPMD> $<$:__TREXIO> $<$:__GREENX> $<$,$>:__FFTW3> diff --git a/src/cp2k_info.F b/src/cp2k_info.F index 8a31af25fe..263c910ff0 100644 --- a/src/cp2k_info.F +++ b/src/cp2k_info.F @@ -274,6 +274,10 @@ CONTAINS flags = TRIM(flags)//" greenx" #endif +#if defined(__OPENPMD) + flags = TRIM(flags)//" openpmd" +#endif + END FUNCTION cp2k_flags ! ************************************************************************************************** diff --git a/src/cp_realspace_grid_openpmd.F b/src/cp_realspace_grid_openpmd.F new file mode 100644 index 0000000000..767b609d85 --- /dev/null +++ b/src/cp_realspace_grid_openpmd.F @@ -0,0 +1,90 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief A wrapper around pw_to_openpmd() which accepts particle_list_type +!> \author Ole Schuett, Franz Poeschel +! ************************************************************************************************** +MODULE cp_realspace_grid_openpmd + USE atomic_kind_types, ONLY: get_atomic_kind + USE kinds, ONLY: dp + USE particle_list_types, ONLY: particle_list_type + USE pw_types, ONLY: pw_r3d_rs_type + USE realspace_grid_openpmd, ONLY: pw_to_openpmd +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + PUBLIC :: cp_pw_to_openpmd + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_realspace_grid_openpmd' + +CONTAINS + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param unit_nr ... +!> \param title ... +!> \param particles ... +!> \param zeff ... +!> \param stride ... +!> \param zero_tails ... +!> \param silent ... +!> \param mpi_io ... +! ************************************************************************************************** + SUBROUTINE cp_pw_to_openpmd( & + pw, & + unit_nr, & + title, & + particles, & + zeff, & + stride, & + zero_tails, & + silent, & + mpi_io & + ) + TYPE(pw_r3d_rs_type), INTENT(IN) :: pw + INTEGER, INTENT(IN) :: unit_nr + CHARACTER(*), INTENT(IN), OPTIONAL :: title + TYPE(particle_list_type), POINTER :: particles + REAL(KIND=dp), DIMENSION(:), OPTIONAL :: zeff + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: stride + LOGICAL, INTENT(IN), OPTIONAL :: zero_tails, silent, mpi_io + + INTEGER :: i, n + INTEGER, ALLOCATABLE, DIMENSION(:) :: particles_z + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: particles_r + TYPE(particle_list_type), POINTER :: my_particles + + NULLIFY (my_particles) + my_particles => particles + IF (ASSOCIATED(my_particles)) THEN + n = my_particles%n_els + ALLOCATE (particles_z(n)) + ALLOCATE (particles_r(3, n)) + DO i = 1, n + CALL get_atomic_kind(my_particles%els(i)%atomic_kind, z=particles_z(i)) + particles_r(:, i) = my_particles%els(i)%r(:) + END DO + + CALL pw_to_openpmd(pw=pw, unit_nr=unit_nr, title=title, & + particles_z=particles_z, particles_r=particles_r, & + particles_zeff=zeff, & + stride=stride, zero_tails=zero_tails, & + silent=silent, mpi_io=mpi_io) + ELSE + CALL pw_to_openpmd(pw=pw, unit_nr=unit_nr, title=title, & + stride=stride, zero_tails=zero_tails, & + silent=silent, mpi_io=mpi_io) + END IF + + END SUBROUTINE cp_pw_to_openpmd + +END MODULE cp_realspace_grid_openpmd diff --git a/src/f77_interface.F b/src/f77_interface.F index 68b7d3423d..36e91e127d 100644 --- a/src/f77_interface.F +++ b/src/f77_interface.F @@ -38,6 +38,7 @@ MODULE f77_interface cp_get_default_logger, cp_logger_create, cp_logger_release, cp_logger_retain, & cp_logger_type, cp_rm_default_logger, cp_to_string USE cp_output_handling, ONLY: cp_iterate + USE cp_output_handling_openpmd, ONLY: cp_openpmd_output_finalize USE cp_result_methods, ONLY: get_results,& test_for_result USE cp_result_types, ONLY: cp_result_type @@ -358,6 +359,7 @@ CONTAINS CALL rm_timer_env() CALL rm_mp_perf_env() CALL string_table_deallocate(0) + CALL cp_openpmd_output_finalize() IF (finalize_mpi) THEN CALL mp_world_finalize() END IF diff --git a/src/input/cp_output_handling_openpmd.F b/src/input/cp_output_handling_openpmd.F new file mode 100644 index 0000000000..674e48cc2b --- /dev/null +++ b/src/input/cp_output_handling_openpmd.F @@ -0,0 +1,740 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief routines to handle the output, The idea is to remove the +!> decision of wheter to output and what to output from the code +!> that does the output, and centralize it here. +!> \note +!> These were originally together with the log handling routines, +!> but have been spawned off. Some dependencies are still there, +!> and some of the comments about log handling also applies to output +!> handling: @see cp_log_handling +! ************************************************************************************************** +MODULE cp_output_handling_openpmd + USE cp_output_handling, only: cp_print_key_should_output, cp_p_file + USE cp_files, ONLY: close_file, & + open_file + USE cp_iter_types, ONLY: cp_iteration_info_release, & + cp_iteration_info_retain, & + cp_iteration_info_type, & + each_desc_labels, & + each_possible_labels + USE cp_log_handling, ONLY: cp_logger_generate_filename, & + cp_logger_get_default_unit_nr, & + cp_logger_get_unit_nr, & + cp_logger_type, & + cp_to_string + USE input_keyword_types, ONLY: keyword_create, & + keyword_release, & + keyword_type + USE input_section_types, ONLY: section_add_keyword, & + section_add_subsection, & + section_create, & + section_release, & + section_type, & + section_vals_get_subs_vals, & + section_vals_type, & + section_vals_val_get + USE kinds, ONLY: default_path_length, & + default_string_length + USE machine, ONLY: m_mov + USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_file_delete, & + mp_file_get_amode, & + mp_file_type +#ifdef __OPENPMD + USE openpmd_api, ONLY: & + openpmd_access_create, & + openpmd_attributable_type, openpmd_iteration_type, openpmd_mesh_type, & + openpmd_particle_species_type, & + openpmd_record_type, & + openpmd_series_create, openpmd_series_type, & + openpmd_type_int, openpmd_json_merge, openpmd_get_default_extension +#endif + USE string_utilities, ONLY: compress, & + s2a +#include "../base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + + LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_output_handling_openpmd' + PUBLIC :: cp_openpmd_print_key_unit_nr, cp_openpmd_print_key_finished_output + PUBLIC :: cp_openpmd_get_value_unit_nr, cp_openpmd_per_call_value_type + PUBLIC :: cp_openpmd_output_finalize + PUBLIC :: cp_openpmd_get_default_extension + PUBLIC :: cp_openpmd_close_iterations + +#ifdef __OPENPMD + TYPE :: cp_openpmd_per_call_value_type + TYPE(openpmd_series_type) :: series + TYPE(openpmd_iteration_type) :: iteration + ! TYPE(openpmd_mesh_type) :: mesh + ! TYPE(openpmd_particle_species_type) :: particle_species + CHARACTER(len=default_string_length) :: name_prefix ! e.g. 'WFN_00008_1' + END TYPE cp_openpmd_per_call_value_type + + TYPE :: cp_openpmd_per_call_type + INTEGER :: key ! unit_nr + TYPE(cp_openpmd_per_call_value_type) :: value + END TYPE cp_openpmd_per_call_type + + TYPE :: cp_current_iteration_counter_type + INTEGER :: flat_iteration = 0 + INTEGER, ALLOCATABLE :: complex_iteration(:) + INTEGER :: complex_iteration_depth = 0 + END TYPE cp_current_iteration_counter_type + + TYPE :: cp_openpmd_per_callsite_value_type + ! openPMD output Series. + TYPE(openpmd_series_type) :: output_series + ! Information on the last Iteration that was written to, including + ! CP2Ks complex Iteration number and its associated contiguous scalar + ! openPMD Iteration number. + TYPE(cp_current_iteration_counter_type) :: iteration_counter + END TYPE cp_openpmd_per_callsite_value_type + + TYPE :: cp_openpmd_per_callsite_type + CHARACTER(len=default_string_length) :: key ! openpmd_basename + TYPE(cp_openpmd_per_callsite_value_type) :: value + END TYPE cp_openpmd_per_callsite_type + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Begin data members for openPMD output. ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Map that associates opened unit numbers with their associated openPMD content. + ! Since CP2K logically opens a new file for every single dataset, multiple + ! unit numbers may point to the same openPMD Series. + TYPE(cp_openpmd_per_call_type), ALLOCATABLE :: cp_openpmd_per_call(:) + INTEGER :: cp_num_openpmd_per_call = 0 + INTEGER :: cp_capacity_openpmd_per_call = 0 + + ! Map that associates callsites from which functions of this module may be invoked + ! to their associated openPMD content. + ! This stores the actual output Series (which stays open across calls from the + ! same callsite) and the Iteration counter (which associates complex CP2k + ! Iterations with flattened scalar Iteration indexes in the openPMD output). + TYPE(cp_openpmd_per_callsite_type), ALLOCATABLE, TARGET :: cp_openpmd_per_callsite(:) + INTEGER :: cp_num_openpmd_per_callsite = 0 + INTEGER :: cp_capacity_openpmd_per_callsite = 0 + + ! This is currently hardcoded, reallocation in case of greater needed map sizes + ! is not (yet) supported. However, the maps should normally not grow to large + ! sizes: + ! + ! * cp_openpmd_per_call will normally contain one single element, since a + ! (virtual) file is opened, written and then closed. + ! The output routines normally do not contain interleaved open-write-close + ! logic. + ! * cp_openpmd_per_callsite will normally contain a handful of elements, + ! equal to the number of output modules activated in the input file + ! (and in openPMD: equal to the number of output Series). + ! There are not 100 of them. + INTEGER, PARAMETER :: cp_allocation_size = 100 + ! Some default settings. May be overwritten / extended by specifying a JSON/TOML + ! config in the input file. + CHARACTER(len=*), PARAMETER :: cp_default_backend_config = & + "[hdf5]"//new_line('a')// & + "# will be overridden by particle flushes"//new_line('a')// & + "independent_stores = false"//new_line('a')// & + "dont_warn_unused_keys = ['independent_stores']"//new_line('a')// & + ""//new_line('a')// & + "[adios2]"//new_line('a')// & + "# discard any attributes written on ranks other than 0"//new_line('a')// & + "attribute_writing_ranks = 0"//new_line('a')// & + "[adios2.engine]"//new_line('a')// & + "# CP2K generally has many small IO operations, "//new_line('a')// & + "# so stage IO memory to the buffer first and then "//new_line('a')// & + "# run it all at once, instead of writing to disk directly."//new_line('a')// & + "# Save memory by specifying 'disk' here instead."//new_line('a')// & + "# TODO: In future, maybe implement some input variable"//new_line('a')// & + "# to specify intervals at which to flush to disk."//new_line('a')// & + "preferred_flush_target = 'buffer'"//new_line('a') +#ifndef _WIN32 + CHARACTER(len=*), PARAMETER :: cp_default_backend_config_non_windows = & + "# Raise the BufferChunkSize to the maximum (2GB), since large operations"//new_line('a')// & + "# improve IO performance and the allocation overhead only cuts into"//new_line('a')// & + "# virtual memory (except on Windows, hence do not do that there)"//new_line('a')// & + "[adios2.engine.parameters]"//new_line('a')// & + "BufferChunkSize = 2147381248"//new_line('a') +#endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! End data members for openPMD output. ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#else ! defined(__OPENPMD) + + TYPE :: cp_openpmd_per_call_value_type + ! nothing there + END TYPE cp_openpmd_per_call_value_type + +#endif + +CONTAINS + +#ifdef __OPENPMD + ! Helper functions for interacting with the two maps declared above. + #:set name_suffixes = ['unit_nr', 'filedata'] + #:set key_types = ['INTEGER', 'CHARACTER(len=default_string_length)'] + #:set value_types = ['cp_openpmd_per_call_value_type', 'cp_openpmd_per_callsite_value_type'] + #:set map_storages = ['cp_openpmd_per_call', 'cp_openpmd_per_callsite'] + #:set map_counters = ['cp_num_openpmd_per_call', 'cp_num_openpmd_per_callsite'] + #:set map_capacities = ['cp_capacity_openpmd_per_call', 'cp_capacity_openpmd_per_callsite'] + + #:for name_suffix, key_type, value_type, map_storage, map_counter, map_capacity in zip(name_suffixes, key_types, value_types, map_storages, map_counters, map_capacities) + + ! TODO No reallocation support for now, change cp_allocation_size if larger sizes needed. + + ! ************************************************************************************************** + !> \brief ... + !> \param key ... + !> \param value ... + ! ************************************************************************************************** + FUNCTION cp_openpmd_add_${name_suffix}$ (key, value) RESULT(index) + ${key_type}$, INTENT(in) :: key + TYPE(${value_type}$), INTENT(in) :: value + INTEGER :: index + + LOGICAL :: check_capacity + INTEGER :: i + + ! Check if the key already exists + DO i = 1, ${map_counter}$ + IF (${map_storage}$ (i)%key == key) THEN + ${map_storage}$ (i)%value = value + index = i + RETURN + END IF + END DO + + IF (${map_capacity}$ == 0) THEN + ALLOCATE (${map_storage}$ (cp_allocation_size)) + ${map_capacity}$ = cp_allocation_size + END IF + + ! No idea how to do reallocations, so for now just assert that they're not needed + check_capacity = ${map_counter}$ < ${map_capacity}$ + CPASSERT(check_capacity) + + ! Add a new entry + ${map_counter}$ = ${map_counter}$+1 + ${map_storage}$ (${map_counter}$)%key = key + ${map_storage}$ (${map_counter}$)%value = value + index = ${map_counter}$ + END FUNCTION cp_openpmd_add_${name_suffix}$ + + ! ************************************************************************************************** + !> \brief ... + !> \param key ... + !> \return ... + ! ************************************************************************************************** + FUNCTION cp_openpmd_get_index_${name_suffix}$ (key) RESULT(index) + ${key_type}$, INTENT(in) :: key + INTEGER :: index + + INTEGER :: i + + index = -1 + + DO i = 1, ${map_counter}$ + IF (${map_storage}$ (i)%key == key) THEN + index = i + RETURN + END IF + END DO + END FUNCTION cp_openpmd_get_index_${name_suffix}$ + + FUNCTION cp_openpmd_get_value_${name_suffix}$ (key) RESULT(value) + ${key_type}$, INTENT(in) :: key + TYPE(${value_type}$) :: value + + INTEGER :: i + + i = cp_openpmd_get_index_${name_suffix}$ (key) + IF (i == -1) RETURN + + value = ${map_storage}$ (i)%value + END FUNCTION cp_openpmd_get_value_${name_suffix}$ + + ! ************************************************************************************************** + !> \brief ... + !> \param key ... + !> \return ... + ! ************************************************************************************************** + FUNCTION cp_openpmd_remove_${name_suffix}$ (key) RESULT(was_found) + ${key_type}$, INTENT(in) :: key + LOGICAL :: was_found + + INTEGER :: i + + was_found = .FALSE. + + DO i = 1, ${map_counter}$ + IF (${map_storage}$ (i)%key == key) THEN + was_found = .TRUE. + IF (i /= ${map_counter}$) THEN + ! Swap last element to now freed place + ${map_storage}$ (i) = ${map_storage}$ (${map_counter}$) + END IF + + ${map_counter}$ = ${map_counter}$-1 + IF (${map_counter}$ == 0) THEN + DEALLOCATE (${map_storage}$) + ${map_capacity}$ = 0 + END IF + RETURN + END IF + END DO + END FUNCTION cp_openpmd_remove_${name_suffix}$ + + #:endfor + +! ************************************************************************************************** +!> \brief Simplified version of cp_print_key_generate_filename. Since an openPMD Series encompasses +! multiple datasets that would be separate outputs in e.g. .cube files, this needs not +! consider dataset names for creation of a filename. +!> \param logger ... +!> \param print_key ... +!> \param openpmd_basename ... +!> \param extension ... +!> \param my_local ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_print_key_generate_openpmd_filename(logger, print_key, openpmd_basename, extension) RESULT(filename) + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), POINTER :: print_key + CHARACTER(len=*), INTENT(IN) :: openpmd_basename, extension + CHARACTER(len=default_path_length) :: filename + + CHARACTER(len=default_path_length) :: outPath, root + CHARACTER(len=default_string_length) :: outName + INTEGER :: my_ind1, my_ind2 + LOGICAL :: has_root + + CALL section_vals_val_get(print_key, "FILENAME", c_val=outPath) + IF (outPath(1:1) == '=') THEN + CPASSERT(LEN(outPath) - 1 <= LEN(filename)) + filename = outPath(2:) + RETURN + END IF + IF (outPath == "__STD_OUT__") outPath = "" + outName = outPath + has_root = .FALSE. + my_ind1 = INDEX(outPath, "/") + my_ind2 = LEN_TRIM(outPath) + IF (my_ind1 /= 0) THEN + has_root = .TRUE. + DO WHILE (INDEX(outPath(my_ind1 + 1:my_ind2), "/") /= 0) + my_ind1 = INDEX(outPath(my_ind1 + 1:my_ind2), "/") + my_ind1 + END DO + IF (my_ind1 == my_ind2) THEN + outName = "" + ELSE + outName = outPath(my_ind1 + 1:my_ind2) + END IF + END IF + + IF (.NOT. has_root) THEN + root = TRIM(logger%iter_info%project_name) + ELSE IF (outName == "") THEN + root = outPath(1:my_ind1)//TRIM(logger%iter_info%project_name) + ELSE + root = outPath(1:my_ind1) + END IF + + filename = ADJUSTL(TRIM(root)//"_"//TRIM(openpmd_basename)//TRIM(extension)) + + END FUNCTION cp_print_key_generate_openpmd_filename + +! ************************************************************************************************** +!> \brief CP2K Iteration numbers are n-dimensional while openPMD Iteration numbers are scalars. +! This checks if the Iteration number has changed from the previous call (stored in +! openpmd_file%iteration_counter) and updates it if needed. +!> \param logger ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_advance_iteration_number(logger, openpmd_file) RESULT(did_advance_iteration) + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_openpmd_per_callsite_value_type) :: openpmd_file + LOGICAL :: did_advance_iteration + + INTEGER :: len + + #:set ic = 'openpmd_file%iteration_counter' + + did_advance_iteration = .FALSE. + len = SIZE(logger%iter_info%iteration) + IF (len /= ${ic}$%complex_iteration_depth) THEN + did_advance_iteration = .TRUE. + ${ic}$%complex_iteration_depth = len + ALLOCATE (${ic}$%complex_iteration(len)) + ELSE + did_advance_iteration & + = ANY(${ic}$%complex_iteration(1:len) & + /= logger%iter_info%iteration(1:len)) + END IF + + IF (.NOT. did_advance_iteration) RETURN + + ${ic}$%flat_iteration = ${ic}$%flat_iteration + 1 + ${ic}$%complex_iteration(1:len) & + = logger%iter_info%iteration(1:len) + + END FUNCTION cp_advance_iteration_number + +! ************************************************************************************************** +!> \brief CP2K deals with output handles in terms of unit numbers. +! The openPMD output logic does not change this association. +! For this, we need to emulate unit numbers as (1) they are not native to openPMD and +! (2) a single openPMD Series might contain multiple datasets treated logically by CP2K +! as distinct outputs. As a result, a single unit number is resolved by the openPMD logic +! to the values represented by the cp_openpmd_per_call_value_type struct, +! containing the output Series and the referred datasets therein (Iteration number, +! name prefix for meshes and particles, referred output Series ...). +!> \param series ... +!> \param middle_name ... +!> \param logger ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_openpmd_create_unit_nr_entry(openpmd_file_index, middle_name, logger) RESULT(res) + INTEGER :: openpmd_file_index + CHARACTER(len=*), INTENT(IN) :: middle_name + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_openpmd_per_call_value_type) :: res + + LOGICAL, SAVE :: opened_new_iteration = .FALSE. + TYPE(openpmd_attributable_type) :: attr + TYPE(cp_openpmd_per_callsite_value_type), POINTER :: opmd + + opmd => cp_openpmd_per_callsite(openpmd_file_index)%value + + res%series = opmd%output_series + + opened_new_iteration = cp_advance_iteration_number(logger, opmd) + + res%iteration = opmd%output_series%write_iteration(opmd%iteration_counter%flat_iteration) + res%name_prefix = TRIM(middle_name) + + IF (opened_new_iteration) THEN + attr = res%iteration%as_attributable() + CALL attr%set_attribute_vec_int( & + "ndim_iteration_index", & + opmd%iteration_counter%complex_iteration) + END IF + END FUNCTION cp_openpmd_create_unit_nr_entry + +! ************************************************************************************************** +!> \brief Check if there is already an output Series created for the callsite identified +! by openpmd_basename. If so, then return it (by index), otherwise open the Series now +! and return the index then. + FUNCTION cp_openpmd_get_openpmd_file_entry(openpmd_basename, filename, openpmd_config, logger, use_mpi) RESULT(file_index) + CHARACTER(len=*), INTENT(IN) :: openpmd_basename, filename, openpmd_config + TYPE(cp_logger_type), POINTER :: logger + LOGICAL :: use_mpi + INTEGER :: file_index + CHARACTER(:), ALLOCATABLE :: merged_config + + CHARACTER(len=default_string_length), SAVE :: basename_copied = ' ' + TYPE(cp_openpmd_per_callsite_value_type) :: emplace_new + + INTEGER :: handle + TYPE(cp_openpmd_per_callsite_value_type) :: series_data + TYPE(openpmd_iteration_type) :: iteration + INTEGER :: i + + basename_copied(1:LEN_TRIM(openpmd_basename)) = TRIM(openpmd_basename) + + file_index = cp_openpmd_get_index_filedata(basename_copied) + + CALL timeset('openpmd_close_iterations', handle) + DO i = 1, cp_num_openpmd_per_callsite + IF (i /= file_index) THEN + series_data = cp_openpmd_per_callsite(i)%value + iteration = series_data%output_series%get_iteration( & + series_data%iteration_counter%flat_iteration) + IF (.NOT. iteration%closed()) THEN + CALL iteration%close() + END IF + END IF + END DO + CALL timestop(handle) + + IF (file_index /= -1) RETURN + +#ifndef _WIN32 + merged_config = openpmd_json_merge(cp_default_backend_config, cp_default_backend_config_non_windows) +#else + merged_config = cp_default_backend_config +#endif + IF (use_mpi) THEN + merged_config = openpmd_json_merge(merged_config, openpmd_config, logger%para_env) + emplace_new%output_series = openpmd_series_create( & + filename, openpmd_access_create, logger%para_env, merged_config) + ELSE + merged_config = openpmd_json_merge(merged_config, openpmd_config) + emplace_new%output_series = openpmd_series_create( & + filename, openpmd_access_create, config=merged_config) + END IF + DEALLOCATE (merged_config) + file_index = cp_openpmd_add_filedata(basename_copied, emplace_new) + END FUNCTION cp_openpmd_get_openpmd_file_entry + +#else ! defined(__OPENPMD) + + FUNCTION cp_openpmd_get_value_unit_nr(key) RESULT(value) + INTEGER, INTENT(in) :: key + TYPE(cp_openpmd_per_call_value_type) :: value + + MARK_USED(key) + MARK_USED(value) + CPABORT("CP2K compiled without the openPMD-api") + + END FUNCTION cp_openpmd_get_value_unit_nr + +#endif + +! ************************************************************************************************** +!> \brief Close all outputs. +! ************************************************************************************************** + SUBROUTINE cp_openpmd_output_finalize() +#ifdef __OPENPMD + INTEGER :: i + DO i = 1, cp_num_openpmd_per_callsite + DEALLOCATE (cp_openpmd_per_callsite(i)%value%iteration_counter%complex_iteration) + CALL cp_openpmd_per_callsite(i)%value%output_series%close() + END DO + IF (ALLOCATED(cp_openpmd_per_callsite)) THEN + DEALLOCATE (cp_openpmd_per_callsite) + END IF + cp_num_openpmd_per_callsite = 0 +#endif + END SUBROUTINE cp_openpmd_output_finalize + +! ************************************************************************************************** +!> \brief ... +!> \param logger ... +!> \param basis_section ... +!> \param print_key_path ... +!> \param extension ... +!> \param middle_name ... +!> \param local ... +!> \param log_filename ... +!> \param ignore_should_output ... +!> \param do_backup ... +!> \param is_new_file true if this rank created a new (or rewound) file, false otherwise +!> \param mpi_io True if the file should be opened in parallel on all processors belonging to +!> the communicator group. Automatically disabled if the file form or access mode +!> is unsuitable for MPI IO. Return value indicates whether MPI was actually used +!> and therefore the flag must also be passed to the file closing directive. +!> \param fout Name of the actual file where the output will be written. Needed mainly for MPI IO +!> because inquiring the filename from the MPI filehandle does not work across +!> all MPI libraries. +!> \param openpmd_basename Used to associate an identifier to each callsite of this module +!> \param use_openpmd ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_openpmd_print_key_unit_nr(logger, basis_section, print_key_path, & + middle_name, ignore_should_output, & + mpi_io, & + fout, openpmd_basename) RESULT(res) + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), INTENT(IN) :: basis_section + CHARACTER(len=*), INTENT(IN), OPTIONAL :: print_key_path + CHARACTER(len=*), INTENT(IN), OPTIONAL :: middle_name + LOGICAL, INTENT(IN), OPTIONAL :: ignore_should_output + LOGICAL, INTENT(INOUT), OPTIONAL :: mpi_io + CHARACTER(len=default_path_length), INTENT(OUT), & + OPTIONAL :: fout + CHARACTER(len=*), INTENT(IN), OPTIONAL :: openpmd_basename + INTEGER :: res + +#ifdef __OPENPMD + + CHARACTER(len=default_path_length) :: filename + + CHARACTER(len=default_string_length) :: openpmd_config, outPath, file_extension + LOGICAL :: found, & + my_mpi_io, & + my_should_output, & + replace + INTEGER :: openpmd_file_index, openpmd_call_index + TYPE(section_vals_type), POINTER :: print_key + + my_mpi_io = .FALSE. + replace = .FALSE. + found = .FALSE. + res = -1 + IF (PRESENT(mpi_io)) THEN +#if defined(__parallel) + IF (logger%para_env%num_pe > 1 .AND. mpi_io) THEN + my_mpi_io = .TRUE. + ELSE + my_mpi_io = .FALSE. + END IF +#else + my_mpi_io = .FALSE. +#endif + ! Set return value + mpi_io = my_mpi_io + END IF + NULLIFY (print_key) + CPASSERT(ASSOCIATED(logger)) + CPASSERT(basis_section%ref_count > 0) + CPASSERT(logger%ref_count > 0) + my_should_output = BTEST(cp_print_key_should_output(logger%iter_info, & + basis_section, print_key_path, used_print_key=print_key), cp_p_file) + IF (PRESENT(ignore_should_output)) my_should_output = my_should_output .OR. ignore_should_output + IF (.NOT. my_should_output) RETURN + IF (logger%para_env%is_source() .OR. my_mpi_io) THEN + + CALL section_vals_val_get(print_key, "FILENAME", c_val=outPath) + CALL section_vals_val_get(print_key, "OPENPMD_EXTENSION", c_val=file_extension) + CALL section_vals_val_get(print_key, "OPENPMD_CFG_FILE", c_val=openpmd_config) + IF (LEN_TRIM(openpmd_config) == 0) THEN + CALL section_vals_val_get(print_key, "OPENPMD_CFG", c_val=openpmd_config) + ELSE + openpmd_config = "@"//openpmd_config + END IF + filename = cp_print_key_generate_openpmd_filename(logger, print_key, openpmd_basename, file_extension) + + IF (PRESENT(fout)) THEN + fout = filename + END IF + + openpmd_file_index = cp_openpmd_get_openpmd_file_entry( & + openpmd_basename, filename, openpmd_config, logger, my_mpi_io) + + OPEN (newunit=res, status='scratch', action='write') + openpmd_call_index = cp_openpmd_add_unit_nr( & + res, & + cp_openpmd_create_unit_nr_entry( & + openpmd_file_index, middle_name, logger)) + + ELSE + res = -1 + END IF +#else + MARK_USED(logger) + MARK_USED(basis_section) + MARK_USED(print_key_path) + MARK_USED(middle_name) + MARK_USED(ignore_should_output) + MARK_USED(mpi_io) + MARK_USED(fout) + MARK_USED(openpmd_basename) + res = 0 + CPABORT("CP2K compiled without the openPMD-api") +#endif + END FUNCTION cp_openpmd_print_key_unit_nr + +! ************************************************************************************************** +!> \brief should be called after you finish working with a unit obtained with +!> cp_openpmd_print_key_unit_nr, so that the file that might have been opened +!> can be closed. +!> +!> the inputs should be exactly the same of the corresponding +!> cp_openpmd_print_key_unit_nr +!> \param unit_nr ... +!> \param logger ... +!> \param basis_section ... +!> \param print_key_path ... +!> \param local ... +!> \param ignore_should_output ... +!> \param mpi_io True if file was opened in parallel with MPI +!> \param use_openpmd ... +!> \note +!> closes if the corresponding filename of the printkey is +!> not __STD_OUT__ +! ************************************************************************************************** + SUBROUTINE cp_openpmd_print_key_finished_output(unit_nr, logger, basis_section, & + print_key_path, local, ignore_should_output, & + mpi_io) + INTEGER, INTENT(INOUT) :: unit_nr + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), INTENT(IN) :: basis_section + CHARACTER(len=*), INTENT(IN), OPTIONAL :: print_key_path + LOGICAL, INTENT(IN), OPTIONAL :: local, ignore_should_output, & + mpi_io + +#ifdef __OPENPMD + + CHARACTER(len=default_string_length) :: outPath + LOGICAL :: my_local, my_mpi_io, & + my_should_output + TYPE(section_vals_type), POINTER :: print_key + + my_local = .FALSE. + my_mpi_io = .FALSE. + NULLIFY (print_key) + IF (PRESENT(local)) my_local = local + IF (PRESENT(mpi_io)) my_mpi_io = mpi_io + CPASSERT(ASSOCIATED(logger)) + CPASSERT(basis_section%ref_count > 0) + CPASSERT(logger%ref_count > 0) + my_should_output = BTEST(cp_print_key_should_output(logger%iter_info, basis_section, & + print_key_path, used_print_key=print_key), cp_p_file) + IF (PRESENT(ignore_should_output)) my_should_output = my_should_output .OR. ignore_should_output + IF (my_should_output .AND. (my_local .OR. & + logger%para_env%is_source() .OR. & + my_mpi_io)) THEN + CALL section_vals_val_get(print_key, "FILENAME", c_val=outPath) + IF (cp_openpmd_remove_unit_nr(unit_nr)) THEN + CLOSE (unit_nr) + END IF + + unit_nr = -1 + END IF + CPASSERT(unit_nr == -1) + unit_nr = -1 +#else + MARK_USED(unit_nr) + MARK_USED(logger) + MARK_USED(basis_section) + MARK_USED(print_key_path) + MARK_USED(local) + MARK_USED(ignore_should_output) + MARK_USED(mpi_io) + CPABORT("CP2K compiled without the openPMD-api") +#endif + END SUBROUTINE cp_openpmd_print_key_finished_output + + SUBROUTINE cp_openpmd_close_iterations() +#ifdef __OPENPMD + INTEGER :: handle + TYPE(cp_openpmd_per_callsite_value_type) :: series_data + TYPE(openpmd_iteration_type) :: iteration + INTEGER :: i + + CALL timeset('openpmd_close_iterations', handle) + DO i = 1, cp_num_openpmd_per_callsite + series_data = cp_openpmd_per_callsite(i)%value + iteration = series_data%output_series%get_iteration( & + series_data%iteration_counter%flat_iteration) + IF (.NOT. iteration%closed()) THEN + CALL iteration%close() + END IF + END DO + CALL timestop(handle) +#endif + END SUBROUTINE cp_openpmd_close_iterations + + FUNCTION cp_openpmd_get_default_extension() RESULT(extension) + CHARACTER(len=default_string_length) :: extension + +#ifdef __OPENPMD + extension = openpmd_get_default_extension() +#else + extension = ".bp5" +#endif + + END FUNCTION cp_openpmd_get_default_extension + +END MODULE cp_output_handling_openpmd diff --git a/src/input_cp2k_print_dft.F b/src/input_cp2k_print_dft.F index d303580108..036f2cc551 100644 --- a/src/input_cp2k_print_dft.F +++ b/src/input_cp2k_print_dft.F @@ -164,6 +164,9 @@ MODULE input_cp2k_print_dft USE qs_mom_types, ONLY: create_mom_section USE string_utilities, ONLY: newline, & s2a + + USE cp_output_handling_openpmd, ONLY: cp_openpmd_get_default_extension + #include "./base/base_uses.f90" IMPLICIT NONE @@ -615,7 +618,23 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) - CALL create_mo_cubes_section(print_key) + CALL create_mo_section(print_key, "MO_CUBES", "cube", [2, 2, 2], "STRIDE 1 1 1", high_print_level, "write_cube") + CALL keyword_create(keyword, __LOCATION__, name="APPEND", & + description="append the cube files when they already exist", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="max_file_size_mb", & + description="Limits the size of the cube file by choosing a suitable stride. Zero means no limit.", & + usage="MAX_FILE_SIZE_MB 1.5", default_r_val=0.0_dp) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + + CALL create_mo_section( & + print_key, "MO_OPENPMD", "openPMD", [1, 1, 1], "STRIDE 2 2 2", debug_print_level + 1, "write_openpmd") + CALL add_generic_openpmd_arguments(print_key) CALL section_add_subsection(section, print_key) CALL section_release(print_key) @@ -700,43 +719,18 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) - CALL cp_print_key_section_create(print_key, __LOCATION__, name="E_DENSITY_CUBE", & - description="Controls the printing of cube files with "// & - "the electronic density and, for LSD calculations, the spin density.", & - print_level=high_print_level, filename="") - CALL keyword_create(keyword, __LOCATION__, name="stride", & - description="The stride (X,Y,Z) used to write the cube file "// & - "(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// & - " 1 number valid for all components.", & - usage="STRIDE 2 2 2", n_var=-1, default_i_vals=[2, 2, 2], type_of_var=integer_t) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="DENSITY_INCLUDE", & - description="Which parts of the density to include. In GAPW the electronic density "// & - "is divided into a hard and a soft component, and the default (TOTAL_HARD_APPROX) "// & - "is to approximate the hard density as a spherical gaussian and to print the smooth "// & - "density accurately. This avoids potential artefacts originating from the hard density. "// & - "If the TOTAL_DENSITY keyword is used the hard density will be computed more accurately "// & - "but may introduce non-physical features. The SOFT_DENSITY keyword will lead to only the "// & - "soft density being printed. In GPW these options have no effect and the cube file will "// & - "only contain the valence electron density.", & - usage="DENSITY_INCLUDE TOTAL_HARD_APPROX", & - enum_c_vals=s2a("TOTAL_HARD_APPROX", "TOTAL_DENSITY", "SOFT_DENSITY"), & - enum_desc=s2a("Print (hard+soft) density where the hard components shape is approximated", & - "Print (hard+soft) density. Only has an effect "// & - "if PAW atoms are present. NOTE: The total "// & - "in real space might exhibit unphysical features "// & - "like spikes due to the finite and thus "// & - "truncated g vector", & - "Print only the soft density"), & - enum_i_vals=[e_dens_total_hard_approx, & - e_dens_total_density, & - e_dens_soft_density], & - default_i_val=e_dens_total_hard_approx) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) + CALL create_e_density_section( & + print_key, & + "E_DENSITY_OPENPMD", & + "openPMD", & + [1, 1, 1], & + "STRIDE 1 1 1", & + debug_print_level + 1) + CALL add_generic_openpmd_arguments(print_key) + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + CALL create_e_density_section(print_key, "E_DENSITY_CUBE", "cube", [2, 2, 2], "STRIDE 2 2 2", high_print_level) CALL keyword_create(keyword, __LOCATION__, name="APPEND", & description="append the cube files when they already exist", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) @@ -952,35 +946,22 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) - CALL cp_print_key_section_create(print_key, __LOCATION__, "ELF_CUBE", & - description="Controls printing of cube files with"// & - " the electron localization function (ELF). Note that"// & - " the value of ELF is defined between 0 and 1: Pauli kinetic energy density normalized"// & - " by the kinetic energy density of a uniform el. gas of same density.", & - print_level=high_print_level, filename="") - CALL keyword_create(keyword, __LOCATION__, name="stride", & - description="The stride (X,Y,Z) used to write the cube file "// & - "(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// & - " 1 number valid for all components.", & - usage="STRIDE 2 2 2", n_var=-1, default_i_vals=[2, 2, 2], type_of_var=integer_t) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) + CALL create_elf_print_section(print_key, "ELF_CUBE", & + "cube", & + [2, 2, 2], "STRIDE 2 2 2", high_print_level, "") CALL keyword_create(keyword, __LOCATION__, name="APPEND", & description="append the cube files when they already exist", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="density_cutoff", & - description=" ", & - usage="density_cutoff 0.0001", & - repeats=.FALSE., & - n_var=1, & - type_of_var=real_t, & - default_r_val=1.0e-10_dp) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + CALL create_elf_print_section(print_key, "ELF_OPENPMD", & + "openPMD", & + [1, 1, 1], "STRIDE 1 1 1", debug_print_level + 1, "") + CALL add_generic_openpmd_arguments(print_key) CALL section_add_subsection(section, print_key) CALL section_release(print_key) @@ -2193,40 +2174,77 @@ CONTAINS END SUBROUTINE create_bandstructure_section + SUBROUTINE add_generic_openpmd_arguments(print_key) + TYPE(section_type), POINTER :: print_key + + TYPE(keyword_type), POINTER :: keyword + NULLIFY (keyword) + + CALL keyword_create(keyword, __LOCATION__, name="OPENPMD_EXTENSION", & + description="Filename extension for openPMD files, including the dot and "// & + "(for optionally activating file encoding) a file expansion pattern.", & + default_c_val="_%06T."//cp_openpmd_get_default_extension(), & + type_of_var=char_t) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="OPENPMD_CFG", & + description="Inline runtime config for openPMD output. Note that inline "// & + "specifications are subject to restrictions imposed by the input "// & + "file format, making this option useful only for very simple use cases. "// & + "Refer to OPENPMD_CFG_FILE for anything else.", & + default_c_val="{}", type_of_var=char_t) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="OPENPMD_CFG_FILE", & + description="Runtime config file for openPMD output. This parameter takes precedence over OPENPMD_CFG.", default_c_val="", & + type_of_var=char_t) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE add_generic_openpmd_arguments + ! ************************************************************************************************** !> \brief creates the input section for dealing with homo lumos, including dumping cubes !> \param print_key ... ! ************************************************************************************************** - SUBROUTINE create_mo_cubes_section(print_key) - TYPE(section_type), POINTER :: print_key + SUBROUTINE create_mo_section( & + print_key, section_name, description, stride_default, stride_usage, & + print_level, do_write_keyname) + + TYPE(section_type), POINTER :: print_key + CHARACTER(len=*), INTENT(IN) :: section_name, description, stride_usage, do_write_keyname + INTEGER, DIMENSION(3), INTENT(IN) :: stride_default + INTEGER, INTENT(IN) :: print_level TYPE(keyword_type), POINTER :: keyword NULLIFY (keyword) - CALL cp_print_key_section_create(print_key, __LOCATION__, "MO_CUBES", & - description="Controls the printing of the molecular orbitals (MOs) as cube files."// & + CALL cp_print_key_section_create(print_key, __LOCATION__, section_name, & + description="Controls the printing of the molecular orbitals (MOs) as " & + //TRIM(ADJUSTL(description))// & + " files."// & " It can be used during a Real Time calculation to print the MOs."// & " In this case, the density corresponding to the time dependent MO is printed"// & " instead of the wave-function.", & - print_level=high_print_level, filename="") + print_level=print_level, filename="") CALL keyword_create(keyword, __LOCATION__, name="stride", & - description="The stride (X,Y,Z) used to write the cube file "// & - "(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// & + description="The stride (X,Y,Z) used to write the "//TRIM(ADJUSTL(description))//" file "// & + "(larger values result in smaller "// & + TRIM(ADJUSTL(description))// & + " files). You can provide 3 numbers (for X,Y,Z) or"// & " 1 number valid for all components.", & - usage="STRIDE 2 2 2", n_var=-1, default_i_vals=[2, 2, 2], type_of_var=integer_t) + usage=stride_usage, n_var=-1, default_i_vals=stride_default, type_of_var=integer_t) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="max_file_size_mb", & - description="Limits the size of the cube file by choosing a suitable stride. Zero means no limit.", & - usage="MAX_FILE_SIZE_MB 1.5", default_r_val=0.0_dp) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="write_cube", & - description="If the MO cube file should be written. If false, the eigenvalues are still computed."// & + CALL keyword_create(keyword, __LOCATION__, name=do_write_keyname, & + description="If the MO " & + //TRIM(ADJUSTL(description)) & + //" file should be written. If false, the eigenvalues are still computed."// & " Can also be useful in combination with STM calculations", & default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(print_key, keyword) @@ -2234,14 +2252,16 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="nlumo", & description="If the printkey is activated controls the number of lumos"// & - " that are printed and dumped as a cube (-1=all)", & + " that are printed and dumped as "//TRIM(ADJUSTL(description))//" (-1=all)", & default_i_val=0) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) CALL keyword_create( & keyword, __LOCATION__, name="nhomo", & - description="If the printkey is activated controls the number of homos that dumped as a cube (-1=all),"// & + description="If the printkey is activated controls the number of homos that dumped as "// & + TRIM(ADJUSTL(description))// & + " (-1=all),"// & " eigenvalues are always all dumped", & default_i_val=1) CALL section_add_keyword(print_key, keyword) @@ -2249,20 +2269,109 @@ CONTAINS CALL keyword_create( & keyword, __LOCATION__, name="homo_list", & - description="If the printkey is activated controls the index of homos dumped as a cube,"// & + description="If the printkey is activated controls the index of homos dumped as openPMD,"// & " eigenvalues are always all dumped. It overrides nhomo.", & usage="HOMO_LIST {integer} {integer} .. {integer} ", type_of_var=integer_t, & n_var=-1, repeats=.TRUE.) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="APPEND", & - description="append the cube files when they already exist", & - default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + END SUBROUTINE create_mo_section + + SUBROUTINE create_e_density_section( & + print_key, section_name, description, stride_default, & + stride_usage, print_level) + + TYPE(section_type), POINTER :: print_key + CHARACTER(len=*), INTENT(IN) :: section_name, description, stride_usage + INTEGER, DIMENSION(3), INTENT(IN) :: stride_default + INTEGER, INTENT(IN) :: print_level + + TYPE(keyword_type), POINTER :: keyword + + NULLIFY (keyword) + + CALL cp_print_key_section_create(print_key, __LOCATION__, name=section_name, & + description="Controls the printing of "//TRIM(ADJUSTL(description))//" files with "// & + "the electronic density and, for LSD calculations, the spin density.", & + print_level=print_level, filename="") + CALL keyword_create(keyword, __LOCATION__, name="stride", & + description="The stride (X,Y,Z) used to write the "//TRIM(ADJUSTL(description))//" file "// & + "(larger values result in smaller "// & + TRIM(ADJUSTL(description))// & + " files). You can provide 3 numbers (for X,Y,Z) or"// & + " 1 number valid for all components.", & + usage=stride_usage, n_var=-1, default_i_vals=stride_default, type_of_var=integer_t) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) - END SUBROUTINE create_mo_cubes_section + CALL keyword_create(keyword, __LOCATION__, name="DENSITY_INCLUDE", & + description="Which parts of the density to include. In GAPW the electronic density "// & + "is divided into a hard and a soft component, and the default (TOTAL_HARD_APPROX) "// & + "is to approximate the hard density as a spherical gaussian and to print the smooth "// & + "density accurately. This avoids potential artefacts originating from the hard density. "// & + "If the TOTAL_DENSITY keyword is used the hard density will be computed more accurately "// & + "but may introduce non-physical features. The SOFT_DENSITY keyword will lead to only the "// & + "soft density being printed. In GPW these options have no effect and the cube file will "// & + "only contain the valence electron density.", & + usage="DENSITY_INCLUDE TOTAL_HARD_APPROX", & + enum_c_vals=s2a("TOTAL_HARD_APPROX", "TOTAL_DENSITY", "SOFT_DENSITY"), & + enum_desc=s2a("Print (hard+soft) density where the hard components shape is approximated", & + "Print (hard+soft) density. Only has an effect "// & + "if PAW atoms are present. NOTE: The total "// & + "in real space might exhibit unphysical features "// & + "like spikes due to the finite and thus "// & + "truncated g vector", & + "Print only the soft density"), & + enum_i_vals=[e_dens_total_hard_approx, & + e_dens_total_density, & + e_dens_soft_density], & + default_i_val=e_dens_total_hard_approx) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_e_density_section + + ! ************************************************************************************************** +!> \brief Helper to create ELF print sections (cube or openPMD) +!> \param print_key section handle (output) +!> \param section_name name of the section (e.g. "ELF_CUBE" or "ELF_OPENPMD") +!> \param description Either "cube" or "openPMD", for the descriptions. +!> \param stride_default default stride values +!> \param stride_usage usage string for stride +!> \param print_level print level +!> \param filename output filename (empty string for default) + SUBROUTINE create_elf_print_section( & + print_key, section_name, description, stride_default, stride_usage, print_level, filename) + + TYPE(section_type), POINTER :: print_key + CHARACTER(len=*), INTENT(IN) :: section_name, description, stride_usage, filename + INTEGER, DIMENSION(3), INTENT(IN) :: stride_default + INTEGER, INTENT(IN) :: print_level + TYPE(keyword_type), POINTER :: keyword + + NULLIFY (keyword) + + CALL cp_print_key_section_create(print_key, __LOCATION__, section_name, & + description="Controls printing of "//TRIM(ADJUSTL(description))// & + " files with the electron localization function (ELF). "// & + "Note that the value of ELF is defined between 0 and 1: "// & + "Pauli kinetic energy density normalized by the kinetic energy density "// & + "of a uniform el. gas of same density.", print_level=print_level, filename=filename) + + CALL keyword_create(keyword, __LOCATION__, name="stride", & + description="The stride (X,Y,Z) used to write the file (larger values result in smaller files). "// & + "You can provide 3 numbers (for X,Y,Z) or 1 number valid for all components.", & + usage=stride_usage, n_var=-1, default_i_vals=stride_default, type_of_var=integer_t) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="density_cutoff", & + description=" ", usage="density_cutoff 0.0001", repeats=.FALSE., n_var=1, & + type_of_var=real_t, default_r_val=1.0e-10_dp) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_elf_print_section ! ************************************************************************************************** !> \brief ... diff --git a/src/mpiwrap/message_passing.F b/src/mpiwrap/message_passing.F index c0f553d504..e08e3d2c92 100644 --- a/src/mpiwrap/message_passing.F +++ b/src/mpiwrap/message_passing.F @@ -1094,13 +1094,14 @@ CONTAINS #endif !$OMP MASTER -#if defined(__DLAF) - ! DLA-Future requires that the MPI library supports THREAD_MULTIPLE mode +#if defined(__DLAF) || defined(__OPENPMD) + ! Both DLA-Future and (some IO backends of) the openPMD-api require + ! that the MPI library supports THREAD_MULTIPLE mode CALL mpi_init_thread(MPI_THREAD_MULTIPLE, provided_tsl, ierr) IF (ierr /= 0) CALL mp_stop(ierr, "mpi_init_thread @ mp_world_init") IF (provided_tsl < MPI_THREAD_MULTIPLE) THEN - CALL mp_stop(0, "MPI library does not support the requested level of threading (MPI_THREAD_MULTIPLE), "// & - "required by DLA-Future. Build CP2K without DLA-Future.") + CALL mp_stop(0, "MPI library does not support the requested level of threading (MPI_THREAD_MULTIPLE),"// & + " required by DLA-Future/openPMD-api. Build CP2K without DLA-Future and openPMD-api.") END IF #else CALL mpi_init_thread(MPI_THREAD_SERIALIZED, provided_tsl, ierr) diff --git a/src/openPMD.F b/src/openPMD.F new file mode 100644 index 0000000000..f8ca81001d --- /dev/null +++ b/src/openPMD.F @@ -0,0 +1,1044 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +MODULE openPMD_api + +#ifdef __OPENPMD + + USE ISO_C_BINDING, ONLY: & + C_ASSOCIATED, & + C_CHAR, & + C_DOUBLE, & + C_F_POINTER, & + C_INT, & + C_INT64_T, & + C_LOC, & + C_NULL_CHAR, & + C_NULL_PTR, & + C_PTR, & + C_SIZE_T + USE kinds, ONLY: default_string_length, dp, sp + USE message_passing, ONLY: mp_comm_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + INTEGER, PARAMETER :: openpmd_access_create = 0 + INTEGER, PARAMETER :: openpmd_access_read_only = 1 + + INTEGER, PARAMETER :: openpmd_type_char = 0 + INTEGER, PARAMETER :: openpmd_type_uchar = 1 + INTEGER, PARAMETER :: openpmd_type_schar = 2 + INTEGER, PARAMETER :: openpmd_type_short = 3 + INTEGER, PARAMETER :: openpmd_type_int = 4 + INTEGER, PARAMETER :: openpmd_type_long = 5 + INTEGER, PARAMETER :: openpmd_type_longlong = 6 + INTEGER, PARAMETER :: openpmd_type_ushort = 7 + INTEGER, PARAMETER :: openpmd_type_uint = 8 + INTEGER, PARAMETER :: openpmd_type_ulong = 9 + INTEGER, PARAMETER :: openpmd_type_ulonglong = 10 + INTEGER, PARAMETER :: openpmd_type_float = 11 + INTEGER, PARAMETER :: openpmd_type_double = 12 + INTEGER, PARAMETER :: openpmd_type_long_double = 13 + INTEGER, PARAMETER :: openpmd_type_cfloat = 14 + INTEGER, PARAMETER :: openpmd_type_cdouble = 15 + INTEGER, PARAMETER :: openpmd_type_clong_double = 16 + INTEGER, PARAMETER :: openpmd_type_string = 17 + INTEGER, PARAMETER :: openpmd_type_vec_char = 18 + INTEGER, PARAMETER :: openpmd_type_vec_short = 19 + INTEGER, PARAMETER :: openpmd_type_vec_int = 20 + INTEGER, PARAMETER :: openpmd_type_vec_long = 21 + INTEGER, PARAMETER :: openpmd_type_vec_longlong = 22 + INTEGER, PARAMETER :: openpmd_type_vec_uchar = 23 + INTEGER, PARAMETER :: openpmd_type_vec_ushort = 24 + INTEGER, PARAMETER :: openpmd_type_vec_uint = 25 + INTEGER, PARAMETER :: openpmd_type_vec_ulong = 26 + INTEGER, PARAMETER :: openpmd_type_vec_ulonglong = 27 + INTEGER, PARAMETER :: openpmd_type_vec_float = 28 + INTEGER, PARAMETER :: openpmd_type_vec_double = 29 + INTEGER, PARAMETER :: openpmd_type_vec_long_double = 30 + INTEGER, PARAMETER :: openpmd_type_vec_cfloat = 31 + INTEGER, PARAMETER :: openpmd_type_vec_cdouble = 32 + INTEGER, PARAMETER :: openpmd_type_vec_clong_double = 33 + INTEGER, PARAMETER :: openpmd_type_vec_schar = 34 + INTEGER, PARAMETER :: openpmd_type_vec_string = 35 + INTEGER, PARAMETER :: openpmd_type_arr_dbl_7 = 36 + INTEGER, PARAMETER :: openpmd_type_bool = 37 + + #:set dataset_types_f = ['REAL(kind=dp)'] + #:set dataset_types_enum = ['double'] + #:set dimensions = [1, 2, 3] + #:set upcasts = [("Series", "Attributable"), ("Iteration", "Attributable"), ("RecordComponent", "Attributable"), ("Mesh", "RecordComponent"), ("Record", "RecordComponent"), ("Mesh", "MeshRecordComponent"), ("MeshRecordComponent", "RecordComponent")] + #:set fortran_names = {"RecordComponent": "record_component", "MeshRecordComponent": "mesh_record_component"} + + TYPE openpmd_attributable_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: set_attribute_vec_int => openpmd_attributable_set_attribute_vec_int + PROCEDURE, PUBLIC :: series_flush => openpmd_attributable_series_flush + END TYPE openpmd_attributable_type + + TYPE openpmd_series_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: as_pointer + PROCEDURE, PUBLIC :: as_attributable => openpmd_series_as_attributable + PROCEDURE, PUBLIC :: close => openpmd_series_close + PROCEDURE, PUBLIC :: present => openpmd_series_present + PROCEDURE, PUBLIC :: write_iteration => openpmd_series_write_iteration + PROCEDURE, PUBLIC :: get_iteration => openpmd_series_get_iteration + END TYPE openpmd_series_type + + TYPE openpmd_iteration_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: as_attributable => openpmd_iteration_as_attributable + PROCEDURE, PUBLIC :: get_mesh => openpmd_iteration_get_mesh + PROCEDURE, PUBLIC :: get_particle_species => openpmd_iteration_get_particle_species + PROCEDURE, PUBLIC :: close => openpmd_iteration_close + PROCEDURE, PUBLIC :: closed => openpmd_iteration_closed + END TYPE openpmd_iteration_type + + TYPE openpmd_mesh_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: as_record_component => openpmd_mesh_as_record_component + PROCEDURE, PUBLIC :: set_axis_labels => openpmd_mesh_set_axis_labels + PROCEDURE, PUBLIC :: set_grid_global_offset => openpmd_mesh_set_grid_global_offset + PROCEDURE, PUBLIC :: set_grid_spacing => openpmd_mesh_set_grid_spacing + PROCEDURE, PUBLIC :: set_position => openpmd_mesh_set_position + END TYPE openpmd_mesh_type + + TYPE openpmd_particle_species_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: get_record => openpmd_particle_species_get_record + END TYPE openpmd_particle_species_type + + TYPE openpmd_record_component_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: make_empty => openpmd_record_component_make_empty + PROCEDURE, PUBLIC :: make_constant_zero => openpmd_record_component_make_constant_zero + PROCEDURE, PUBLIC :: reset_dataset => openpmd_record_component_reset_dataset + #:for type_f, type_enum in zip(dataset_types_f, dataset_types_enum) + #:for dim in dimensions + PROCEDURE, PRIVATE :: store_chunk_${dim}$d_${type_enum}$ & + => openpmd_record_component_store_chunk_${dim}$d_${type_enum}$ + PROCEDURE, PUBLIC :: store_chunk_span_${dim}$d_${type_enum}$ & + => openpmd_record_component_store_chunk_span_${dim}$d_${type_enum}$ + #:endfor + #:endfor + GENERIC, PUBLIC :: store_chunk => & + #:set zipped = zip(dataset_types_f, dataset_types_enum) + #:set n = len([x for x in zipped]) + #:set m = len(dimensions) + #:set count_up_to = m * n + #:for i, type_f, type_enum in zip(range(n), dataset_types_f, dataset_types_enum) + #:for j, dim in enumerate(dimensions) + #:if i * n + j == count_up_to - 1 + store_chunk_${dim}$d_${type_enum}$ + #:else + store_chunk_${dim}$d_${type_enum}$, & + #:endif + #:endfor + #:endfor + #:del count_up_to + #:del zipped + #:del n + #:del m + END TYPE openpmd_record_component_type + + TYPE openpmd_mesh_record_component_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + END TYPE openpmd_mesh_record_component_type + + TYPE openpmd_record_type + PRIVATE + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + CONTAINS + PROCEDURE, PUBLIC :: as_record_component => openpmd_record_as_record_component + PROCEDURE, PUBLIC :: get_component => openpmd_record_get_component + END TYPE openpmd_record_type + + #:for dim in dimensions + TYPE openpmd_dynamic_memory_view_type_${dim}$d + PRIVATE + INTEGER, DIMENSION(${dim}$) :: chunk_extent + TYPE(C_PTR) :: c_ptr = C_NULL_PTR + + CONTAINS + + #:for type_f, type_enum in zip(dataset_types_f, dataset_types_enum) + PROCEDURE, PUBLIC :: resolve_${type_enum}$ & + => openpmd_dynamic_memory_view_resolve_${dim}$d_${type_enum}$ + #:endfor + END TYPE openpmd_dynamic_memory_view_type_${dim}$d + #:endfor + + ! Types + PUBLIC :: openpmd_attributable_type + PUBLIC :: openpmd_iteration_type + PUBLIC :: openpmd_mesh_type + PUBLIC :: openpmd_particle_species_type + PUBLIC :: openpmd_record_component_type + PUBLIC :: openpmd_record_type + PUBLIC :: openpmd_series_type, openpmd_series_create + + ! Helpers + PUBLIC :: openpmd_access_create, openpmd_access_read_only + PUBLIC :: openpmd_get_default_extension + PUBLIC :: openpmd_json_merge + + ! openPMD datatypes + PUBLIC :: openpmd_type_arr_dbl_7 + PUBLIC :: openpmd_type_bool + PUBLIC :: openpmd_type_cdouble + PUBLIC :: openpmd_type_cfloat + PUBLIC :: openpmd_type_char + PUBLIC :: openpmd_type_clong_double + PUBLIC :: openpmd_type_double + PUBLIC :: openpmd_type_float + PUBLIC :: openpmd_type_int + PUBLIC :: openpmd_type_long + PUBLIC :: openpmd_type_long_double + PUBLIC :: openpmd_type_longlong + PUBLIC :: openpmd_type_schar + PUBLIC :: openpmd_type_short + PUBLIC :: openpmd_type_string + PUBLIC :: openpmd_type_uchar + PUBLIC :: openpmd_type_uint + PUBLIC :: openpmd_type_ulong + PUBLIC :: openpmd_type_ulonglong + PUBLIC :: openpmd_type_ushort + PUBLIC :: openpmd_type_vec_cdouble + PUBLIC :: openpmd_type_vec_cfloat + PUBLIC :: openpmd_type_vec_char + PUBLIC :: openpmd_type_vec_clong_double + PUBLIC :: openpmd_type_vec_double + PUBLIC :: openpmd_type_vec_float + PUBLIC :: openpmd_type_vec_int + PUBLIC :: openpmd_type_vec_long + PUBLIC :: openpmd_type_vec_long_double + PUBLIC :: openpmd_type_vec_longlong + PUBLIC :: openpmd_type_vec_schar + PUBLIC :: openpmd_type_vec_short + PUBLIC :: openpmd_type_vec_string + PUBLIC :: openpmd_type_vec_uchar + PUBLIC :: openpmd_type_vec_uint + PUBLIC :: openpmd_type_vec_ulong + PUBLIC :: openpmd_type_vec_ulonglong + PUBLIC :: openpmd_type_vec_ushort + + #:for child, parent in upcasts + #:set child_f = fortran_names[child] if child in fortran_names else child + #:set parent_f = fortran_names[parent] if parent in fortran_names else parent + PUBLIC :: openpmd_${child_f}$_as_${parent_f}$ + #:endfor + #:for dim in dimensions + PUBLIC :: openpmd_dynamic_memory_view_type_${dim}$d + #:endfor + +#endif + + CONTAINS + +#ifdef __OPENPMD + +! ************************************************************************************************** +!> \brief ... +!> \param this ... +!> \return ... +! ************************************************************************************************** + FUNCTION as_pointer(this) RESULT(res) + CLASS(openpmd_series_type), INTENT(IN) :: this + TYPE(C_PTR) :: res + + res = this%c_ptr + END FUNCTION as_pointer + + #:for child, parent in upcasts + #:set child_f = fortran_names[child] if child in fortran_names else child + #:set parent_f = fortran_names[parent] if parent in fortran_names else parent + FUNCTION openpmd_${child_f}$_as_${parent_f}$ (this) RESULT(res) + CLASS(openpmd_${child_f}$_type), INTENT(IN) :: this + TYPE(openpmd_${parent_f}$_type) :: res + + INTERFACE + SUBROUTINE openpmd_c_${child_f}$_as_${parent_f}$ (child, parent) & + BIND(C, name="openPMD_${child}$_upcast_to_${parent}$") + IMPORT :: C_PTR + TYPE(C_PTR), VALUE :: child + TYPE(C_PTR) :: parent + END SUBROUTINE openpmd_c_${child_f}$_as_${parent_f}$ + END INTERFACE + + CALL openpmd_c_${child_f}$_as_${parent_f}$ (this%c_ptr, res%c_ptr) + END FUNCTION openpmd_${child_f}$_as_${parent_f}$ + #:endfor + +! ************************************************************************************************** +!> \brief ... +!> \param series ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_series_present(this) RESULT(res) + CLASS(openpmd_series_type), INTENT(IN) :: this + LOGICAL :: res + INTEGER :: res_internal + + INTERFACE + FUNCTION openpmd_c_series_present(series) RESULT(res) & + BIND(C, NAME="openPMD_Series_present") + IMPORT :: C_PTR, C_INT + TYPE(C_PTR), VALUE, INTENT(IN) :: series + INTEGER(kind=C_INT) :: res + END FUNCTION openpmd_c_series_present + END INTERFACE + + IF (.NOT. C_ASSOCIATED(this%c_ptr)) THEN + res = .FALSE. + RETURN + END IF + + res_internal = openpmd_c_series_present(this%c_ptr) + res = res_internal /= 0 + END FUNCTION openpmd_series_present + +! ************************************************************************************************** +!> \brief ... +!> \param path ... +!> \param access ... +!> \param mpi_comm ... +!> \param config ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_series_create(path, access, mpi_comm, config) RESULT(series) + CHARACTER(len=*), INTENT(IN) :: path + INTEGER(kind=C_INT), INTENT(IN) :: access + CLASS(mp_comm_type), INTENT(in), OPTIONAL :: mpi_comm + CHARACTER(len=*), INTENT(IN), OPTIONAL, TARGET :: config + TYPE(openpmd_series_type) :: series + + TYPE(c_ptr), SAVE :: comm_c = C_NULL_PTR + CHARACTER(len=default_string_length), TARGET, SAVE :: default_config = "{}" + CHARACTER(:), POINTER :: my_config + INTERFACE +! ************************************************************************************************** +!> \brief ... +!> \param comm_f ... +!> \param comm_c ... +!> \param C ... +!> \param name="CP2K_MPI_Comm_f2c" ... +! ************************************************************************************************** + SUBROUTINE cp2k_c_mpi_comm_f2c(comm_f, comm_c) bind(C, name="CP2K_MPI_Comm_f2c") + import :: c_ptr, mp_comm_type, c_int + INTEGER(kind=C_INT), VALUE :: comm_f + TYPE(c_ptr) :: comm_c + END SUBROUTINE cp2k_c_mpi_comm_f2c + END INTERFACE + INTERFACE + SUBROUTINE openpmd_c_series_create_mpi(path, access, mpi_comm, config, series) & + BIND(C, name="openPMD_Series_create_mpi") + IMPORT :: C_PTR, C_CHAR, C_INT + CHARACTER(kind=C_CHAR), DIMENSION(*) :: path + INTEGER(kind=C_INT), VALUE :: access + TYPE(C_PTR), value :: mpi_comm + CHARACTER(kind=C_CHAR), DIMENSION(*) :: config + TYPE(C_PTR) :: series + END SUBROUTINE openpmd_c_series_create_mpi + END INTERFACE + INTERFACE + SUBROUTINE openpmd_c_series_create(path, access, config, series) & + BIND(C, name="openPMD_Series_create") + IMPORT :: C_PTR, C_CHAR, C_INT + TYPE(C_PTR) :: series + CHARACTER(kind=C_CHAR), DIMENSION(*) :: path + INTEGER(kind=C_INT), VALUE :: access + CHARACTER(kind=C_CHAR), DIMENSION(*) :: config + END SUBROUTINE openpmd_c_series_create + END INTERFACE + + IF (PRESENT(config)) THEN + my_config => config + ELSE + my_config => default_config + END IF + + CPASSERT(.NOT. C_ASSOCIATED(series%c_ptr)) + IF (PRESENT(mpi_comm)) THEN + CALL cp2k_c_mpi_comm_f2c(mpi_comm%get_handle(), comm_c) + CALL openpmd_c_series_create_mpi( & + path=TRIM(path)//C_NULL_CHAR, & + access=access, & + mpi_comm=comm_c, & + config=my_config//C_NULL_CHAR, & + series=series%c_ptr) + ELSE + CALL openpmd_c_series_create( & + path=TRIM(path)//C_NULL_CHAR, & + access=access, & + config=my_config//C_NULL_CHAR, & + series=series%c_ptr) + END IF + CPASSERT(C_ASSOCIATED(series%c_ptr)) + END FUNCTION openpmd_series_create + +! ************************************************************************************************** +!> \brief ... +!> \param series ... +! ************************************************************************************************** + SUBROUTINE openpmd_series_close(series) + CLASS(openpmd_series_type), INTENT(INOUT) :: series + + INTERFACE +! ************************************************************************************************** +!> \brief ... +!> \param series ... +! ************************************************************************************************** + SUBROUTINE openpmd_c_series_close(series) & + BIND(C, name="openPMD_Series_close") + IMPORT :: C_PTR + TYPE(C_PTR), VALUE :: series + END SUBROUTINE openpmd_c_series_close + END INTERFACE + + CPASSERT(C_ASSOCIATED(series%c_ptr)) + CALL openpmd_c_series_close(series=series%c_ptr) + series%c_ptr = C_NULL_PTR + END SUBROUTINE openpmd_series_close + +! ************************************************************************************************** +!> \brief ... +!> \param series ... +!> \param index ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_series_write_iteration(this, index) RESULT(iteration) + CLASS(openpmd_series_type) :: this + INTEGER, TARGET :: index + TYPE(openpmd_iteration_type) :: iteration + + INTEGER(8) :: cast_index + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_series_write_iteration(series, index, iteration) RESULT(status) & + BIND(C, name="openPMD_Series_write_Iteration") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T + TYPE(C_PTR), VALUE :: series + INTEGER(kind=C_INT64_T), VALUE :: index + TYPE(C_PTR) :: iteration + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_series_write_iteration + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(iteration%c_ptr)) + cast_index = index + status = openpmd_c_series_write_iteration(this%c_ptr, cast_index, iteration%c_ptr) + CPASSERT(C_ASSOCIATED(iteration%c_ptr)) + END FUNCTION openpmd_series_write_iteration + + FUNCTION openpmd_series_get_iteration(this, index) RESULT(iteration) + CLASS(openpmd_series_type) :: this + INTEGER, TARGET :: index + TYPE(openpmd_iteration_type) :: iteration + + INTEGER(8) :: cast_index + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_series_get_iteration(series, index, iteration) RESULT(status) & + BIND(C, name="openPMD_Series_get_Iteration") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T + TYPE(C_PTR), VALUE :: series + INTEGER(kind=C_INT64_T), VALUE :: index + TYPE(C_PTR) :: iteration + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_series_get_iteration + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(iteration%c_ptr)) + cast_index = index + status = openpmd_c_series_get_iteration(this%c_ptr, cast_index, iteration%c_ptr) + CPASSERT(C_ASSOCIATED(iteration%c_ptr)) + END FUNCTION openpmd_series_get_iteration + +! ************************************************************************************************** +!> \brief ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_get_default_extension() RESULT(extension) + CHARACTER(len=default_string_length) :: extension + + CHARACTER(kind=C_CHAR), POINTER :: c_extension + INTEGER(C_SIZE_T) :: length_of_c_string, i + INTERFACE + FUNCTION openpmd_c_get_default_extension() RESULT(extension) & + BIND(c, name="openPMD_get_default_extension") + IMPORT :: C_CHAR + CHARACTER(kind=C_CHAR), POINTER :: extension + END FUNCTION openpmd_c_get_default_extension + END INTERFACE + INTERFACE + FUNCTION get_strlen(c_string) RESULT(strlen) & + BIND(C, NAME="CP2K_strlen") + IMPORT :: C_CHAR, C_SIZE_T + CHARACTER(kind=C_CHAR), DIMENSION(*), INTENT(IN) :: c_string + INTEGER(C_SIZE_T) :: strlen + END FUNCTION get_strlen + END INTERFACE + + c_extension => openpmd_c_get_default_extension() + + ! Find the length of the C string + IF (ASSOCIATED(c_extension)) THEN + !c_string_array => c_extension + length_of_c_string = get_strlen(c_extension) + ELSE + RETURN + END IF + + ! Ensure we do not exceed the length of extension + IF (length_of_c_string > default_string_length) THEN + length_of_c_string = default_string_length + END IF + + ! Copy the contents from c_extension to extension + extension = TRIM(ADJUSTL(TRANSFER(c_extension(1:length_of_c_string), extension))) + + DO i = length_of_c_string + 1, LEN(extension) + extension(i:i) = ' ' + END DO + + END FUNCTION openpmd_get_default_extension + +! ************************************************************************************************** +!> \brief ... +!> \param attributable ... +!> \param name ... +!> \param vec ... +! ************************************************************************************************** + SUBROUTINE openpmd_attributable_set_attribute_vec_int(this, name, vec) + CLASS(openpmd_attributable_type) :: this + CHARACTER(len=*), INTENT(IN) :: name + INTEGER, DIMENSION(:), TARGET :: vec + + INTEGER, TARGET :: size_ + INTERFACE + SUBROUTINE openpmd_c_attributable_set_attribute_vec_int(attributable, name, vec, length) & + BIND(C, NAME="openPMD_attributable_set_attribute_vec_int") + IMPORT :: C_CHAR, C_SIZE_T, C_INT, C_PTR + TYPE(C_PTR), VALUE :: attributable + CHARACTER(kind=C_CHAR), DIMENSION(*) :: name + TYPE(C_PTR), VALUE :: vec + INTEGER(kind=C_INT), VALUE :: length + END SUBROUTINE openpmd_c_attributable_set_attribute_vec_int + END INTERFACE + + size_ = SIZE(vec) + CALL openpmd_c_attributable_set_attribute_vec_int( & + this%c_ptr, name//C_NULL_CHAR, C_LOC(vec(1)), size_) + + END SUBROUTINE openpmd_attributable_set_attribute_vec_int + + SUBROUTINE openpmd_attributable_series_flush(this, backendconfig) + CLASS(openpmd_attributable_type) :: this + CHARACTER(len=*), INTENT(IN), OPTIONAL :: backendconfig + CHARACTER(len=default_string_length), SAVE :: my_backendconfig = "" + + INTERFACE + SUBROUTINE openpmd_c_attributable_series_flush(attributable, config) & + BIND(C, NAME="openPMD_Attributable_series_flush") + IMPORT :: C_CHAR, C_SIZE_T, C_INT, C_PTR + TYPE(C_PTR), VALUE :: attributable + CHARACTER(kind=C_CHAR), DIMENSION(*) :: config + END SUBROUTINE openpmd_c_attributable_series_flush + END INTERFACE + + IF (PRESENT(backendconfig)) my_backendconfig = backendconfig + CALL openpmd_c_attributable_series_flush(this%c_ptr, my_backendconfig//C_NULL_CHAR) + END SUBROUTINE openpmd_attributable_series_flush + + FUNCTION openpmd_json_merge(into, from, mpi_comm) RESULT(merged) + CHARACTER(len=*), INTENT(IN) :: into, from + CLASS(mp_comm_type), INTENT(in), OPTIONAL :: mpi_comm + CHARACTER(:), ALLOCATABLE :: merged + CHARACTER(kind=C_CHAR), POINTER :: c_merged + INTEGER(C_SIZE_T) :: length_of_c_string + TYPE(c_ptr), SAVE :: comm_c = C_NULL_PTR + + INTERFACE +! ************************************************************************************************** +!> \brief ... +!> \param comm_f ... +!> \param comm_c ... +!> \param C ... +!> \param name="CP2K_MPI_Comm_f2c" ... +! ************************************************************************************************** + SUBROUTINE cp2k_c_mpi_comm_f2c(comm_f, comm_c) bind(C, name="CP2K_MPI_Comm_f2c") + import :: c_ptr, mp_comm_type, c_int + INTEGER(kind=C_INT), VALUE :: comm_f + TYPE(c_ptr) :: comm_c + END SUBROUTINE cp2k_c_mpi_comm_f2c + END INTERFACE + INTERFACE + FUNCTION openpmd_c_json_merge(into, from, mpi_comm) RESULT(merged) & + BIND(C, NAME="openPMD_json_merge") + IMPORT :: C_CHAR, C_INT, C_PTR + CHARACTER(kind=C_CHAR), DIMENSION(*) :: into, from + TYPE(C_PTR), value :: mpi_comm + CHARACTER(kind=C_CHAR), POINTER :: merged + END FUNCTION openpmd_c_json_merge + END INTERFACE + INTERFACE + FUNCTION get_strlen(c_string) RESULT(strlen) & + BIND(C, NAME="CP2K_strlen") + IMPORT :: C_CHAR, C_SIZE_T + CHARACTER(kind=C_CHAR), DIMENSION(*), INTENT(IN) :: c_string + INTEGER(C_SIZE_T) :: strlen + END FUNCTION get_strlen + END INTERFACE + INTERFACE + SUBROUTINE free(c_string) & + BIND(C, NAME="CP2K_free") + IMPORT :: C_CHAR + CHARACTER(kind=C_CHAR), DIMENSION(*) :: c_string + END SUBROUTINE free + END INTERFACE + + IF (PRESENT(mpi_comm)) THEN + CALL cp2k_c_mpi_comm_f2c(mpi_comm%get_handle(), comm_c) + END IF + c_merged => openpmd_c_json_merge(into//C_NULL_CHAR, from//C_NULL_CHAR, comm_c) + length_of_c_string = get_strlen(c_merged) + ALLOCATE (character(length_of_c_string) :: merged) + merged = TRIM(ADJUSTL(TRANSFER(c_merged(1:length_of_c_string), MERGED))) + CALL free(c_merged) + END FUNCTION openpmd_json_merge + +! ************************************************************************************************** +!> \brief ... +!> \param iteration ... +!> \param name ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_iteration_get_mesh(this, name) RESULT(mesh) + CLASS(openpmd_iteration_type) :: this + CHARACTER(len=*), INTENT(IN) :: name + TYPE(openpmd_mesh_type) :: mesh + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_iteration_get_mesh(iteration, name, mesh) RESULT(status) & + BIND(C, name="openPMD_Iteration_get_mesh") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: iteration + CHARACTER(kind=C_CHAR), DIMENSION(*) :: name + TYPE(C_PTR) :: mesh + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_iteration_get_mesh + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(mesh%c_ptr)) + status = openpmd_c_iteration_get_mesh(this%c_ptr, name//C_NULL_CHAR, mesh%c_ptr) + CPASSERT(C_ASSOCIATED(mesh%c_ptr)) + END FUNCTION openpmd_iteration_get_mesh + +! ************************************************************************************************** +!> \brief ... +!> \param iteration ... +!> \param name ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_iteration_get_particle_species(this, name) RESULT(particle_species) + CLASS(openpmd_iteration_type) :: this + CHARACTER(len=*), INTENT(IN) :: name + TYPE(openpmd_particle_species_type) :: particle_species + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_iteration_get_particle_species(iteration, name, particle_species) RESULT(status) & + BIND(C, name="openPMD_Iteration_get_particle_species") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: iteration + CHARACTER(kind=C_CHAR), DIMENSION(*) :: name + TYPE(C_PTR) :: particle_species + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_iteration_get_particle_species + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(particle_species%c_ptr)) + status = openpmd_c_iteration_get_particle_species( & + this%c_ptr, name//C_NULL_CHAR, particle_species%c_ptr) + CPASSERT(C_ASSOCIATED(particle_species%c_ptr)) + END FUNCTION openpmd_iteration_get_particle_species + +! ************************************************************************************************** +!> \brief ... +!> \param iteration ... +! ************************************************************************************************** + SUBROUTINE openpmd_iteration_close(iteration) + CLASS(openpmd_iteration_type), INTENT(INOUT) :: iteration + + INTERFACE +! ************************************************************************************************** +!> \brief ... +!> \param iteration ... +! ************************************************************************************************** + SUBROUTINE openpmd_c_iteration_close(iteration) & + BIND(C, name="openPMD_Iteration_close") + IMPORT :: C_PTR + TYPE(C_PTR), VALUE :: iteration + END SUBROUTINE openpmd_c_iteration_close + END INTERFACE + + CPASSERT(C_ASSOCIATED(iteration%c_ptr)) + CALL openpmd_c_iteration_close(iteration=iteration%c_ptr) + iteration%c_ptr = C_NULL_PTR + END SUBROUTINE openpmd_iteration_close + + FUNCTION openpmd_iteration_closed(iteration) RESULT(res) + CLASS(openpmd_iteration_type), INTENT(IN) :: iteration + LOGICAL :: res + + INTEGER(kind=C_INT) :: res_tmp + + INTERFACE +! ************************************************************************************************** +!> \brief ... +!> \param iteration ... +! ************************************************************************************************** + FUNCTION openpmd_c_iteration_closed(iteration) RESULT(res) & + BIND(C, name="openPMD_Iteration_closed") + IMPORT :: C_PTR, C_INT + TYPE(C_PTR), VALUE :: iteration + INTEGER(kind=C_INT) :: res + END FUNCTION openpmd_c_iteration_closed + END INTERFACE + + CPASSERT(C_ASSOCIATED(iteration%c_ptr)) + res_tmp = openpmd_c_iteration_closed(iteration=iteration%c_ptr) + IF (res_tmp == 0) THEN + res = .FALSE. + ELSE + res = .TRUE. + END IF + END FUNCTION openpmd_iteration_closed + +! ************************************************************************************************** +!> \brief ... +!> \param record_component ... +!> \param dtype ... +!> \param dimensionality ... +! ************************************************************************************************** + SUBROUTINE openpmd_record_component_make_empty(this, dtype, dimensionality) + CLASS(openpmd_record_component_type) :: this + INTEGER(kind=C_INT) :: dtype + INTEGER :: dimensionality + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_record_component_make_empty(record_component, dtype, dimensionality) RESULT(status) & + BIND(C, name="openPMD_RecordComponent_makeEmpty") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record_component + INTEGER(KIND=C_INT), VALUE :: dtype + INTEGER(KIND=C_INT), VALUE :: dimensionality + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_component_make_empty + END INTERFACE + + status = openpmd_c_record_component_make_empty(this%c_ptr, dtype, dimensionality) + END SUBROUTINE openpmd_record_component_make_empty + +! ************************************************************************************************** +!> \brief ... +!> \param record_component ... +!> \param dtype ... +!> \param extent ... +! ************************************************************************************************** + SUBROUTINE openpmd_record_component_make_constant_zero(this, dtype, extent) + CLASS(openpmd_record_component_type) :: this + INTEGER(kind=C_INT) :: dtype + INTEGER, DIMENSION(:), TARGET :: extent + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_record_component_make_constant( & + record_component, dtype, dimensions, extent, invert, value) RESULT(status) & + BIND(C, name="openPMD_RecordComponent_makeConstant") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record_component + INTEGER(KIND=C_INT), VALUE :: dtype + INTEGER(KIND=C_INT), VALUE :: dimensions, invert + TYPE(C_PTR), VALUE :: extent + TYPE(C_PTR), VALUE :: value + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_component_make_constant + END INTERFACE + + status = openpmd_c_record_component_make_constant( & + this%c_ptr, dtype, SIZE(extent), C_LOC(extent(1)), 1, C_NULL_PTR) + END SUBROUTINE openpmd_record_component_make_constant_zero + + SUBROUTINE openpmd_mesh_set_axis_labels(this, labels) + CLASS(openpmd_mesh_type) :: this + CHARACTER(len=*), DIMENSION(:) :: labels + + INTEGER :: i + CHARACTER(len=default_string_length), ALLOCATABLE, TARGET :: labels_as_null_terminated_strings(:) + TYPE(C_PTR), ALLOCATABLE, TARGET :: labels_as_c_strings(:) + + INTERFACE + SUBROUTINE openpmd_c_mesh_set_axis_labels(mesh, labels, len_labels, invert) & + BIND(C, NAME="openPMD_Mesh_set_axis_labels") + IMPORT :: C_CHAR, C_SIZE_T, C_INT, C_PTR + TYPE(C_PTR), VALUE :: mesh + TYPE(C_PTR), VALUE :: labels + INTEGER(kind=C_INT), VALUE :: len_labels, invert + END SUBROUTINE openpmd_c_mesh_set_axis_labels + END INTERFACE + + ALLOCATE (labels_as_null_terminated_strings(SIZE(labels))) + ALLOCATE (labels_as_c_strings(SIZE(labels))) + + DO i = 1, SIZE(labels) + labels_as_null_terminated_strings(i) = TRIM(labels(i))//C_NULL_CHAR + labels_as_c_strings(i) = C_LOC(labels_as_null_terminated_strings(i)) + END DO + + CALL openpmd_c_mesh_set_axis_labels(this%c_ptr, C_LOC(labels_as_c_strings(1)), SIZE(labels), invert=1) + + DEALLOCATE (labels_as_c_strings) + DEALLOCATE (labels_as_null_terminated_strings) + END SUBROUTINE openpmd_mesh_set_axis_labels + + SUBROUTINE openpmd_record_component_reset_dataset(this, dtype, extent, cfg) + CLASS(openpmd_record_component_type) :: this + INTEGER(kind=C_INT) :: dtype + INTEGER, DIMENSION(:), TARGET :: extent + CHARACTER(len=*), OPTIONAL :: cfg + + CHARACTER(len=default_string_length), SAVE :: my_config = "{}" + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_record_component_reset_dataset( & + record_component, dtype, dimensions, extent, invert, cfg) RESULT(status) & + BIND(C, name="openPMD_RecordComponent_resetDataset") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record_component + INTEGER(KIND=C_INT), VALUE :: dtype + INTEGER(KIND=C_INT), VALUE :: dimensions, invert + TYPE(C_PTR), VALUE :: extent + CHARACTER(kind=C_CHAR), DIMENSION(*) :: cfg + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_component_reset_dataset + END INTERFACE + + IF (PRESENT(cfg)) my_config = cfg + + status = openpmd_c_record_component_reset_dataset( & + this%c_ptr, dtype, SIZE(extent), C_LOC(extent(1)), 1, my_config//C_NULL_CHAR) + END SUBROUTINE openpmd_record_component_reset_dataset + + #:for type_f, type_enum in zip(dataset_types_f, dataset_types_enum) + #:for dim in [1, 2, 3] + ! e.g. ':,:,:' + #:set dimensionality_wildcard = ','.join(':' for _ in range(dim)) + SUBROUTINE openpmd_record_component_store_chunk_${dim}$d_${type_enum}$ (record_component, data_, offset) + CLASS(openpmd_record_component_type) :: record_component + ${type_f}$, TARGET :: data_(${dimensionality_wildcard}$) + INTEGER, DIMENSION(${dim}$), TARGET :: offset, extent + + INTEGER(kind=C_INT) :: status + + INTERFACE + FUNCTION openpmd_c_record_component_store_chunk( & + record_component, dtype, dimensions, offset, extent, invert, data_) RESULT(status) & + BIND(C, name="openPMD_RecordComponent_storeChunk") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record_component + INTEGER(KIND=C_INT), VALUE :: dtype + INTEGER(KIND=C_INT), VALUE :: dimensions, invert + TYPE(C_PTR), VALUE :: offset, extent + TYPE(C_PTR), VALUE :: data_ + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_component_store_chunk + END INTERFACE + + extent = SHAPE(data_) + status = openpmd_c_record_component_store_chunk( & + record_component%c_ptr, & + openpmd_type_${type_enum}$, & + ${dim}$, & + C_LOC(offset(1)), & + C_LOC(extent(1)), & + 1, & + C_LOC(data_)) + + END SUBROUTINE openpmd_record_component_store_chunk_${dim}$d_${type_enum}$ + + FUNCTION openpmd_record_component_store_chunk_span_${dim}$d_${type_enum}$ ( & + record_component, offset, extent) RESULT(memory_view) + CLASS(openpmd_record_component_type) :: record_component + INTEGER, DIMENSION(${dim}$), TARGET :: offset, extent + + TYPE(openpmd_dynamic_memory_view_type_${dim}$d) :: memory_view + + INTEGER(kind=C_INT) :: status + + INTERFACE + FUNCTION openpmd_c_record_component_store_chunk_span( & + record_component, dtype, dimensions, offset, extent, invert, memory_view) RESULT(status) & + BIND(C, name="openPMD_RecordComponent_storeChunkSpan") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record_component + INTEGER(KIND=C_INT), VALUE :: dtype + INTEGER(KIND=C_INT), VALUE :: dimensions, invert + TYPE(C_PTR), VALUE :: offset, extent + TYPE(C_PTR) :: memory_view + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_component_store_chunk_span + END INTERFACE + + memory_view%chunk_extent = extent + status = openpmd_c_record_component_store_chunk_span( & + record_component%c_ptr, & + openpmd_type_${type_enum}$, & + ${dim}$, & + C_LOC(offset(1)), & + C_LOC(extent(1)), & + 1, & + memory_view%c_ptr) + + END FUNCTION openpmd_record_component_store_chunk_span_${dim}$d_${type_enum}$ + + FUNCTION openpmd_dynamic_memory_view_resolve_${dim}$d_${type_enum}$ (memory_view, deallocate) RESULT(buffer) + CLASS(openpmd_dynamic_memory_view_type_${dim}$d) :: memory_view + LOGICAL :: deallocate + ${type_f}$, POINTER :: buffer(${dimensionality_wildcard}$) + + INTEGER(kind=C_INT) :: status, deallocate_c + TYPE(C_PTR), TARGET :: c_buffer + + INTERFACE + FUNCTION openpmd_c_dynamic_memory_view_resolve(memory_view, deallocate, write_buffer) RESULT(status) & + BIND(C, name="openPMD_DynamicMemoryView_resolve") + IMPORT :: C_PTR, C_INT + TYPE(C_PTR), VALUE :: memory_view + INTEGER(KIND=C_INT), VALUE :: deallocate + TYPE(C_PTR) :: write_buffer + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_dynamic_memory_view_resolve + END INTERFACE + + IF (deallocate) THEN + deallocate_c = 1 + ELSE + deallocate_c = 0 + END IF + + status = openpmd_c_dynamic_memory_view_resolve(memory_view%c_ptr, deallocate_c, c_buffer) + CALL c_f_pointer(c_buffer, buffer, memory_view%chunk_extent) + END FUNCTION openpmd_dynamic_memory_view_resolve_${dim}$d_${type_enum}$ + #:endfor + #:endfor + +! ************************************************************************************************** +!> \brief ... +!> \param particle_species ... +!> \param name ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_particle_species_get_record(this, name) RESULT(record) + CLASS(openpmd_particle_species_type) :: this + CHARACTER(len=*), INTENT(IN) :: name + TYPE(openpmd_record_type) :: record + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_particle_species_get_record(particle_species, name, record) RESULT(status) & + BIND(C, name="openPMD_ParticleSpecies_get_Record") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: particle_species + CHARACTER(kind=C_CHAR), DIMENSION(*) :: name + TYPE(C_PTR) :: record + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_particle_species_get_record + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(record%c_ptr)) + status = openpmd_c_particle_species_get_record(this%c_ptr, name//C_NULL_CHAR, record%c_ptr) + CPASSERT(C_ASSOCIATED(record%c_ptr)) + END FUNCTION openpmd_particle_species_get_record + +! ************************************************************************************************** +!> \brief ... +!> \param record ... +!> \param name ... +!> \return ... +! ************************************************************************************************** + FUNCTION openpmd_record_get_component(this, name) RESULT(record_component) + CLASS(openpmd_record_type) :: this + CHARACTER(len=*), INTENT(IN) :: name + TYPE(openpmd_record_component_type) :: record_component + + INTEGER(kind=C_INT) :: status + INTERFACE + FUNCTION openpmd_c_record_get_component(record, name, record_component) RESULT(status) & + BIND(C, name="openPMD_Record_get_Component") + IMPORT :: C_PTR, C_INT64_T, C_INT, C_SIZE_T, C_CHAR + TYPE(C_PTR), VALUE :: record + CHARACTER(kind=C_CHAR), DIMENSION(*) :: name + TYPE(C_PTR) :: record_component + INTEGER(KIND=C_INT) :: status + END FUNCTION openpmd_c_record_get_component + END INTERFACE + + CPASSERT(.NOT. C_ASSOCIATED(record_component%c_ptr)) + status = openpmd_c_record_get_component(this%c_ptr, name//C_NULL_CHAR, record_component%c_ptr) + CPASSERT(C_ASSOCIATED(record_component%c_ptr)) + END FUNCTION openpmd_record_get_component + + #:set f_routines = ["openpmd_mesh_set_grid_global_offset", "openpmd_mesh_set_grid_spacing", "openpmd_mesh_set_position"] + #:set c_functions = ["openPMD_Mesh_setGridGlobalOffset", "openPMD_Mesh_setGridSpacing", "openPMD_Mesh_setPosition"] + #:set f_vector_types = ["REAL(KIND=dp)" for _ in range(3)] + + #:for f_routine, c_function, f_vector_type in zip(f_routines, c_functions, f_vector_types) + SUBROUTINE ${f_routine}$ (this, vec) + CLASS(openpmd_mesh_type) :: this + ${f_vector_type}$, DIMENSION(:), TARGET :: vec + + INTEGER, TARGET :: size_ + INTERFACE + SUBROUTINE openpmd_c_set_vector_attribute(mesh, vec, length, invert) & + BIND(C, NAME="${c_function}$") + IMPORT :: C_CHAR, C_SIZE_T, C_INT, C_PTR + TYPE(C_PTR), VALUE :: mesh + TYPE(C_PTR), VALUE :: vec + INTEGER(kind=C_INT), VALUE :: length, invert + END SUBROUTINE openpmd_c_set_vector_attribute + END INTERFACE + + size_ = SIZE(vec) + CALL openpmd_c_set_vector_attribute(this%c_ptr, C_LOC(vec(1)), size_, 1) + + END SUBROUTINE ${f_routine}$ + #:endfor + +#endif + ! In two lines of this file, `make pretty` applies a wrong indentation and keeps it for the rest of this file. + ! Hence, we end up a bit further to the right than we should. + END MODULE openPMD_api diff --git a/src/openPMD.cpp b/src/openPMD.cpp new file mode 100644 index 0000000000..4c585cadfd --- /dev/null +++ b/src/openPMD.cpp @@ -0,0 +1,1003 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright 2000-2026 CP2K developers group */ +/* */ +/* SPDX-License-Identifier: GPL-2.0-or-later */ +/*----------------------------------------------------------------------------*/ + +// Single-file library, logically containing 4 files: +// openPMD.h openPMD.cpp auxiliary.h auxiliary.c +// Implements C bindings for a subselection of the openPMD-api. + +#ifdef __OPENPMD + +/////////////// +// openPMD.h // +/////////////// + +// #ifndef CP2K_OPENPMD_H +// #define CP2K_OPENPMD_H + +#include + +#include + +#ifdef __cplusplus +extern "C" { +#endif // defined(__cplusplus) + +typedef enum { openPMD_Access_create, openPMD_Access_read_only } openPMD_Access; + +typedef enum { + openPMD_Type_CHAR, + openPMD_Type_UCHAR, + openPMD_Type_SCHAR, + openPMD_Type_SHORT, + openPMD_Type_INT, + openPMD_Type_LONG, + openPMD_Type_LONGLONG, + openPMD_Type_USHORT, + openPMD_Type_UINT, + openPMD_Type_ULONG, + openPMD_Type_ULONGLONG, + openPMD_Type_FLOAT, + openPMD_Type_DOUBLE, + openPMD_Type_LONG_DOUBLE, + openPMD_Type_CFLOAT, + openPMD_Type_CDOUBLE, + openPMD_Type_CLONG_DOUBLE, + openPMD_Type_STRING, + openPMD_Type_VEC_CHAR, + openPMD_Type_VEC_SHORT, + openPMD_Type_VEC_INT, + openPMD_Type_VEC_LONG, + openPMD_Type_VEC_LONGLONG, + openPMD_Type_VEC_UCHAR, + openPMD_Type_VEC_USHORT, + openPMD_Type_VEC_UINT, + openPMD_Type_VEC_ULONG, + openPMD_Type_VEC_ULONGLONG, + openPMD_Type_VEC_FLOAT, + openPMD_Type_VEC_DOUBLE, + openPMD_Type_VEC_LONG_DOUBLE, + openPMD_Type_VEC_CFLOAT, + openPMD_Type_VEC_CDOUBLE, + openPMD_Type_VEC_CLONG_DOUBLE, + openPMD_Type_VEC_SCHAR, + openPMD_Type_VEC_STRING, + openPMD_Type_ARR_DBL_7, + openPMD_Type_BOOL +} openPMD_Datatype; + +typedef struct openPMD_Attributable_opaque *openPMD_Attributable; + +typedef struct openPMD_Series_opaque *openPMD_Series; + +typedef struct openPMD_Iteration_opaque *openPMD_Iteration; + +typedef struct openPMD_Mesh_opaque *openPMD_Mesh; + +typedef struct openPMD_ParticleSpecies_opaque *openPMD_ParticleSpecies; + +typedef struct openPMD_RecordComponent_opaque *openPMD_MeshRecordComponent; + +typedef struct openPMD_RecordComponent_opaque *openPMD_RecordComponent; + +typedef struct openPMD_Record_opaque *openPMD_Record; + +typedef struct openPMD_DynamicMemoryView_opaque *openPMD_DynamicMemoryView; + +// Actually uint64_t, but ISO C bindings for Fortran only have that as a GNU +// extension, so we use signed integers and convert +typedef uint64_t openPMD_Iteration_Index_t; + +/******************* + * Series members. * + *******************/ + +int openPMD_Series_create( + // in + char const *filename, openPMD_Access access, char const *config, + // out + openPMD_Series *series); + +int openPMD_Series_create_mpi( + // in + char const *filename, openPMD_Access access, MPI_Comm comm, + char const *config, + // out + openPMD_Series *series); + +int openPMD_Series_close( + // in + openPMD_Series series); + +int openPMD_Series_write_Iteration( + // in + openPMD_Series series, openPMD_Iteration_Index_t index, + // out + openPMD_Iteration *iteration); + +int openPMD_Series_get_Iteration( + // in + openPMD_Series series, openPMD_Iteration_Index_t index, + // out + openPMD_Iteration *iteration); + +int openPMD_Series_upcast_to_Attributable( + // in + openPMD_Series series, + // out + openPMD_Attributable *attr); + +int openPMD_Series_present(openPMD_Series series); + +char const *openPMD_get_default_extension(); + +/************************* + * Attributable members. * + *************************/ + +int openPMD_attributable_set_attribute_vec_int(openPMD_Attributable attr, + char const *attr_name, + int const *begin, int length); + +int openPMD_attributable_set_attribute_vec_string(openPMD_Attributable attr, + char const *attr_name, + char const **begin, + int length); + +int openPMD_Attributable_series_flush(openPMD_Attributable attr, char const *); + +/********************** + * Iteration members. * + **********************/ + +int openPMD_Iteration_upcast_to_Attributable( + // in + openPMD_Iteration iteration, + // out + openPMD_Attributable *attr); + +int openPMD_Iteration_get_mesh( + // in + openPMD_Iteration iteration, char const *name, + // out + openPMD_Mesh *mesh); + +int openPMD_Iteration_get_particle_species( + // in + openPMD_Iteration iteration, char const *name, + // out + openPMD_ParticleSpecies *particle_species); + +int openPMD_Iteration_close( + // in + openPMD_Iteration iteration); + +int openPMD_Iteration_closed( + // in + openPMD_Iteration iteration); + +/**************** + * Mesh members * + ****************/ + +int openPMD_Mesh_upcast_to_RecordComponent( + // in + openPMD_Mesh mesh, + // out + openPMD_RecordComponent *rc); + +int openPMD_Mesh_upcast_to_MeshRecordComponent( + // in + openPMD_Mesh mesh, + // out + openPMD_MeshRecordComponent *mrc); + +int openPMD_Mesh_set_axis_labels( + // in + openPMD_Mesh mesh, char const **labels, int len_labels, int invert); + +int openPMD_Mesh_setGridGlobalOffset( + // in + openPMD_Mesh mesh, double const *labels, int len_labels, int invert); + +int openPMD_Mesh_setGridSpacing( + // in + openPMD_Mesh mesh, double const *labels, int len_labels, int invert); + +int openPMD_Mesh_setPosition( + // in + openPMD_Mesh mesh, double const *labels, int len_labels, int invert); + +/*************************** + * RecordComponent members * + ***************************/ + +int openPMD_RecordComponent_upcast_to_Attributable( + // in + openPMD_RecordComponent rc, + // out + openPMD_Attributable *attr); + +int openPMD_RecordComponent_resetDataset( + // in + openPMD_RecordComponent rc, openPMD_Datatype, int dimensions, + int const *extent, int invert, char const *cfg); + +int openPMD_RecordComponent_makeEmpty( + // in + openPMD_RecordComponent rc, openPMD_Datatype, int dimensions); + +int openPMD_RecordComponent_makeConstant( + // in + openPMD_RecordComponent rc, openPMD_Datatype, int dimensions, + int const *extent, int invert, void const *value); + +int openPMD_RecordComponent_storeChunk( + // in + openPMD_RecordComponent rc, openPMD_Datatype, int dimensions, + int const *offset, int const *extent, int invert, void *data); + +int openPMD_RecordComponent_storeChunkSpan( + // in + openPMD_RecordComponent rc, openPMD_Datatype, int dimensions, + int const *offset, int const *extent, int invert, + // out + openPMD_DynamicMemoryView *); + +int openPMD_DynamicMemoryView_resolve( + // in + openPMD_DynamicMemoryView, int deallocate, + // out + void **write_buffer); + +/******************************* + * MeshRecordComponent members * + *******************************/ + +int openPMD_MeshRecordComponent_upcast_to_RecordComponent( + // in + openPMD_MeshRecordComponent mrc, + // out + openPMD_RecordComponent *rc); + +/*************************** + * ParticleSpecies members * + ***************************/ + +int openPMD_ParticleSpecies_get_Record( + // in + openPMD_ParticleSpecies, char const *name, + // out + openPMD_Record *); + +/****************** + * Record members * + ******************/ + +int openPMD_Record_upcast_to_RecordComponent( + // in + openPMD_Record record, + // out + openPMD_RecordComponent *rc); + +int openPMD_Record_get_Component( + // in + openPMD_Record record, char const *name, + // out + openPMD_RecordComponent *rc); + +/*********** + * Helpers * + ***********/ + +char *openPMD_json_merge(char const *into, char const *from, + MPI_Comm maybe_comm); + +#ifdef __cplusplus +} // extern "C" +#endif // defined(__cplusplus) + +// #endif // !defined(CP2K_OPENPMD_H) + +///////////////// +// openPMD.cpp // +///////////////// + +#include +#if !OPENPMDAPI_VERSION_GE(0, 17, 0) +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace implementation { +namespace { +struct DynamicMemoryView { + openPMD_Datatype dtype; + std::any memory_view; +}; + +constexpr auto datatype_c_to_cxx(openPMD_Datatype dt) -> openPMD::Datatype { + switch (dt) { + case openPMD_Type_CHAR: + return openPMD::Datatype::CHAR; + case openPMD_Type_UCHAR: + return openPMD::Datatype::UCHAR; + case openPMD_Type_SCHAR: + return openPMD::Datatype::SCHAR; + case openPMD_Type_SHORT: + return openPMD::Datatype::SHORT; + case openPMD_Type_INT: + return openPMD::Datatype::INT; + case openPMD_Type_LONG: + return openPMD::Datatype::LONG; + case openPMD_Type_LONGLONG: + return openPMD::Datatype::LONGLONG; + case openPMD_Type_USHORT: + return openPMD::Datatype::USHORT; + case openPMD_Type_UINT: + return openPMD::Datatype::UINT; + case openPMD_Type_ULONG: + return openPMD::Datatype::ULONG; + case openPMD_Type_ULONGLONG: + return openPMD::Datatype::ULONGLONG; + case openPMD_Type_FLOAT: + return openPMD::Datatype::FLOAT; + case openPMD_Type_DOUBLE: + return openPMD::Datatype::DOUBLE; + case openPMD_Type_LONG_DOUBLE: + return openPMD::Datatype::LONG_DOUBLE; + case openPMD_Type_CFLOAT: + return openPMD::Datatype::CFLOAT; + case openPMD_Type_CDOUBLE: + return openPMD::Datatype::CDOUBLE; + case openPMD_Type_CLONG_DOUBLE: + return openPMD::Datatype::CLONG_DOUBLE; + case openPMD_Type_STRING: + return openPMD::Datatype::STRING; + case openPMD_Type_VEC_CHAR: + return openPMD::Datatype::VEC_CHAR; + case openPMD_Type_VEC_SHORT: + return openPMD::Datatype::VEC_SHORT; + case openPMD_Type_VEC_INT: + return openPMD::Datatype::VEC_INT; + case openPMD_Type_VEC_LONG: + return openPMD::Datatype::VEC_LONG; + case openPMD_Type_VEC_LONGLONG: + return openPMD::Datatype::VEC_LONGLONG; + case openPMD_Type_VEC_UCHAR: + return openPMD::Datatype::VEC_UCHAR; + case openPMD_Type_VEC_USHORT: + return openPMD::Datatype::VEC_USHORT; + case openPMD_Type_VEC_UINT: + return openPMD::Datatype::VEC_UINT; + case openPMD_Type_VEC_ULONG: + return openPMD::Datatype::VEC_ULONG; + case openPMD_Type_VEC_ULONGLONG: + return openPMD::Datatype::VEC_ULONGLONG; + case openPMD_Type_VEC_FLOAT: + return openPMD::Datatype::VEC_FLOAT; + case openPMD_Type_VEC_DOUBLE: + return openPMD::Datatype::VEC_DOUBLE; + case openPMD_Type_VEC_LONG_DOUBLE: + return openPMD::Datatype::VEC_LONG_DOUBLE; + case openPMD_Type_VEC_CFLOAT: + return openPMD::Datatype::VEC_CFLOAT; + case openPMD_Type_VEC_CDOUBLE: + return openPMD::Datatype::VEC_CDOUBLE; + case openPMD_Type_VEC_CLONG_DOUBLE: + return openPMD::Datatype::VEC_CLONG_DOUBLE; + case openPMD_Type_VEC_SCHAR: + return openPMD::Datatype::VEC_SCHAR; + case openPMD_Type_VEC_STRING: + return openPMD::Datatype::VEC_STRING; + case openPMD_Type_ARR_DBL_7: + return openPMD::Datatype::ARR_DBL_7; + case openPMD_Type_BOOL: + return openPMD::Datatype::BOOL; + } + return openPMD::Datatype::UNDEFINED; +} + +constexpr auto access_c_to_cxx(openPMD_Access access) -> openPMD::Access { + switch (access) { + case openPMD_Access_create: + return openPMD::Access::CREATE; + case openPMD_Access_read_only: + return openPMD::Access::READ_ONLY; + } + // unreachable + return static_cast(0); +} + +template +auto Series_create(openPMD_Series *series_param, + Args &&...constructor_args) -> int { + auto series = reinterpret_cast(series_param); + try { + *series = new openPMD::Series(std::forward(constructor_args)...); + } catch (std::exception const &e) { + std::cout << "[Series_create] Caught error: '" << e.what() << "'\n"; + delete *series; + return 1; + } catch (...) { + std::cout << "[Series_create] Caught unknown error.\n"; + delete *series; + return 1; + } + return 0; +} + +/* + * Use to static_cast<> from a C++ subclass From to its parent To. + * FromOpaque and ToOpaque are the corresponding opaque types from + * the C header, pointers need to be reinterpret_cast<>ed to/from the + * actual underlying C++ types. + */ +template +int do_upcast(FromOpaque from_param, ToOpaque *to_param) { + auto from = reinterpret_cast(from_param); + auto to = reinterpret_cast(to_param); + auto from_upcasted = static_cast(from); + *to = from_upcasted; + return 0; +} + +template +auto pointer_to_vector(T *ptr, size_t len, bool invert) { + using resolved_type = std::conditional_t, + std::remove_cv_t, res_t>; + if (!invert) { + return std::vector{ptr, ptr + len}; + } else { + return std::vector{std::reverse_iterator(ptr + len), + std::reverse_iterator(ptr)}; + } +} + +template +auto pointer_to_vector(T *ptr, size_t len, int invert) { + return pointer_to_vector(ptr, len, static_cast(invert)); +} + +#if !OPENPMDAPI_VERSION_GE(0, 17, 0) +auto resolveFilename(char const *unparsed) -> std::optional { + char const *current_char = unparsed; + // First, trim front + while (std::isspace(*current_char) && *current_char != '\0' && + *current_char != '@') { + ++current_char; + } + if (*current_char == '@') { + // Now, trim back + auto begin_filename = current_char + 1; + auto end_filename = begin_filename; + size_t i = 0; + while (*end_filename != '\0') { + ++end_filename; + ++i; + if (i > FILENAME_MAX) { + throw std::runtime_error( + "resolveFilename: Longer than maximum allowed filename."); + } + } + // *end_filename is now equal to '\0' + --end_filename; + while (std::isspace(*end_filename)) { + --end_filename; + } + // need past-the-end iterator + ++end_filename; + return std::string(begin_filename, end_filename); + } else { + return std::nullopt; + } +} + +auto resolveFileContent(char const *maybeStringMaybeFilename, + std::optional maybe_comm) + -> std::variant { + auto filename = resolveFilename(maybeStringMaybeFilename); + if (!filename.has_value()) { + return maybeStringMaybeFilename; + } + auto serialImplementation = [&]() { + std::fstream handle; + handle.open(*filename, std::ios_base::in); + std::stringstream stream; + stream << handle.rdbuf(); + if (!handle.good()) { + throw std::runtime_error("Failed acessing file '" + *filename + "'."); + } + handle.close(); + return stream.str(); + }; +#if openPMD_HAVE_MPI + auto parallelImplementation = [&](MPI_Comm comm) { + return openPMD::auxiliary::collective_file_read(*filename, comm); + }; + if (maybe_comm.has_value()) { + return parallelImplementation(*maybe_comm); + } else { + return serialImplementation(); + } +#else + return serialImplementation(); +#endif +} +#endif + +struct RecordComponent_makeConstant { + template + static int call(openPMD::RecordComponent &rc, void const *value_in) { + if (value_in) { + auto value = static_cast(value_in); + rc.makeConstant(*value); + } else { + rc.makeConstant(Type{}); + } + return 0; + } + static constexpr char const *errorMsg = "RecordComponent_makeConstant"; +}; + +struct RecordComponent_storeChunk { + template + static int call(openPMD::RecordComponent &rc, openPMD::Offset const &o, + openPMD::Extent const &e, void *data) { + rc.storeChunkRaw(static_cast(data), o, e); + return 0; + } + static constexpr char const *errorMsg = "RecordComponent_storeChunk"; +}; + +struct RecordComponent_storeChunkSpan { + template + static int call(openPMD::RecordComponent &rc, openPMD::Offset const &o, + openPMD::Extent const &e, std::any &memory_view) { + auto buffer = rc.storeChunk(o, e); + memory_view = + std::make_any>(std::move(buffer)); + return 0; + } + static constexpr char const *errorMsg = "RecordComponent_storeChunkSpan"; +}; + +struct DynamicMemoryView_resolve { + template static void *call(DynamicMemoryView &memory_view) { + auto &get_buffer = std::any_cast &>( + memory_view.memory_view); + auto span = get_buffer.currentBuffer(); + return span.data(); + } + static constexpr char const *errorMsg = "DynamicMemoryView_resolve"; +}; +} // namespace +} // namespace implementation + +extern "C" { +int openPMD_Series_create(char const *filename, openPMD_Access access, + char const *config, openPMD_Series *series_param) { + if (!config) { + config = ""; + } + return implementation::Series_create( + series_param, filename, implementation::access_c_to_cxx(access), config); +} + +int openPMD_Series_create_mpi(char const *filename, openPMD_Access access, + MPI_Comm comm, char const *config, + openPMD_Series *series_param) { + if (!config) { + config = ""; + } + return implementation::Series_create(series_param, filename, + implementation::access_c_to_cxx(access), + comm, config); +} + +int openPMD_Series_close(openPMD_Series series_param) { + auto series = reinterpret_cast(series_param); + series->close(); + delete series; + return 0; +} + +int openPMD_Series_write_Iteration(openPMD_Series series_param, + openPMD_Iteration_Index_t index, + openPMD_Iteration *iteration) { + auto series = reinterpret_cast(series_param); + auto &res_iteration = series->writeIterations()[index]; + *reinterpret_cast(iteration) = &res_iteration; + return 0; +} + +int openPMD_Series_get_Iteration(openPMD_Series series_param, + openPMD_Iteration_Index_t index, + openPMD_Iteration *iteration) { + auto series = reinterpret_cast(series_param); + auto &res_iteration = series->iterations[index]; + *reinterpret_cast(iteration) = &res_iteration; + return 0; +} + +int openPMD_Series_upcast_to_Attributable( + // in + openPMD_Series series, + // out + openPMD_Attributable *attr) { + return implementation::do_upcast( + series, attr); +} + +int openPMD_Series_present(openPMD_Series series_param) { + auto series = reinterpret_cast(series_param); + bool res = series && *series; + return res ? 1 : 0; +} + +char const *openPMD_get_default_extension() { + constexpr static char const *preferred_order[] = {"bp5", "bp4", "bp", "h5", + "json"}; + auto const &available_extensions = openPMD::getFileExtensions(); + for (auto const &ext : preferred_order) { + if (std::find(available_extensions.begin(), available_extensions.end(), + ext) != available_extensions.end()) { + return ext; + } + } + return nullptr; +} + +int openPMD_attributable_set_attribute_vec_int(openPMD_Attributable attr, + char const *attr_name, + const int *begin, int length) { + auto attributable = reinterpret_cast(attr); + attributable->setAttribute(attr_name, + std::vector{begin, begin + size_t(length)}); + return 0; +} + +int openPMD_attributable_set_attribute_vec_string(openPMD_Attributable attr, + char const *attr_name, + char const **begin, + int length) { + auto attributable = reinterpret_cast(attr); + attributable->setAttribute( + attr_name, std::vector{begin, begin + size_t(length)}); + return 0; +} + +int openPMD_Attributable_series_flush(openPMD_Attributable attr, + char const *config) { + auto attributable = reinterpret_cast(attr); + attributable->seriesFlush(config ? config : ""); + return 0; +} + +int openPMD_Iteration_upcast_to_Attributable( + // in + openPMD_Iteration iteration, + // out + openPMD_Attributable *attr) { + return implementation::do_upcast( + iteration, attr); +} + +int openPMD_Iteration_get_mesh( + // in + openPMD_Iteration iteration_param, char const *name, + // out + openPMD_Mesh *mesh) { + auto iteration = reinterpret_cast(iteration_param); + auto &res = iteration->meshes[name]; + *reinterpret_cast(mesh) = &res; + return 0; +} + +int openPMD_Iteration_get_particle_species( + // in + openPMD_Iteration iteration_param, char const *name, + // out + openPMD_ParticleSpecies *particle_species) { + auto iteration = reinterpret_cast(iteration_param); + auto &res = iteration->particles[name]; + *reinterpret_cast(particle_species) = &res; + return 0; +} + +int openPMD_Iteration_close(openPMD_Iteration iteration_param) { + auto iteration = reinterpret_cast(iteration_param); + iteration->close(); + return 0; +} + +int openPMD_Iteration_closed(openPMD_Iteration iteration_param) { + auto iteration = reinterpret_cast(iteration_param); + return iteration->closed(); +} + +int openPMD_Mesh_upcast_to_RecordComponent( + // in + openPMD_Mesh mesh_param, + // out + openPMD_RecordComponent *rc) { + return implementation::do_upcast( + mesh_param, rc); +} + +int openPMD_Mesh_upcast_to_MeshRecordComponent( + // in + openPMD_Mesh mesh, + // out + openPMD_MeshRecordComponent *mrc) { + return implementation::do_upcast( + mesh, mrc); +} + +int openPMD_Mesh_set_axis_labels( + // in + openPMD_Mesh mesh_param, char const **labels, int len_labels, int invert) { + auto mesh = reinterpret_cast(mesh_param); + mesh->setAxisLabels(implementation::pointer_to_vector( + labels, len_labels, invert)); + return 0; +} + +int openPMD_Mesh_setGridGlobalOffset( + // in + openPMD_Mesh mesh_param, double const *labels, int len_labels, int invert) { + auto mesh = reinterpret_cast(mesh_param); + mesh->setGridGlobalOffset( + implementation::pointer_to_vector(labels, len_labels, invert)); + return 0; +} + +int openPMD_Mesh_setGridSpacing( + // in + openPMD_Mesh mesh_param, double const *labels, int len_labels, int invert) { + auto mesh = reinterpret_cast(mesh_param); + mesh->setGridSpacing( + implementation::pointer_to_vector(labels, len_labels, invert)); + return 0; +} + +int openPMD_Mesh_setPosition( + // in + openPMD_Mesh mesh_param, double const *labels, int len_labels, int invert) { + auto mesh = reinterpret_cast(mesh_param); + mesh->setPosition( + implementation::pointer_to_vector(labels, len_labels, invert)); + return 0; +} + +int openPMD_MeshRecordComponent_upcast_to_RecordComponent( + // in + openPMD_MeshRecordComponent mrc, + // out + openPMD_RecordComponent *rc) { + return implementation::do_upcast(mrc, rc); +} + +int openPMD_RecordComponent_upcast_to_Attributable( + // in + openPMD_RecordComponent rc, + // out + openPMD_Attributable *attr) { + return implementation::do_upcast(rc, attr); +} + +int openPMD_RecordComponent_resetDataset( + // in + openPMD_RecordComponent rc_param, openPMD_Datatype dtype, int dimensions, + int const *extent, int invert, char const *cfg) { + auto rc = reinterpret_cast(rc_param); + openPMD::Dataset ds( + implementation::datatype_c_to_cxx(dtype), + implementation::pointer_to_vector( + extent, dimensions, invert), + cfg ? cfg : "{}"); + rc->resetDataset(std::move(ds)); + return 0; +} + +int openPMD_RecordComponent_makeEmpty( + // in + openPMD_RecordComponent rc_param, openPMD_Datatype dtype, int dimensions) { + auto rc = reinterpret_cast(rc_param); + rc->makeEmpty(implementation::datatype_c_to_cxx(dtype), dimensions); + return 0; +} + +int openPMD_RecordComponent_makeConstant( + // in + openPMD_RecordComponent rc_param, openPMD_Datatype dt, int dimensions, + int const *extent, int invert, void const *value) { + openPMD_RecordComponent_resetDataset(rc_param, dt, dimensions, extent, invert, + NULL); + auto rc = reinterpret_cast(rc_param); + + openPMD::switchDatasetType( + implementation::datatype_c_to_cxx(dt), *rc, value); + + return 0; +} + +int openPMD_RecordComponent_storeChunk( + // in + openPMD_RecordComponent rc_param, openPMD_Datatype dt, int dimensions, + int const *offset, int const *extent, int invert, void *data) { + auto rc = reinterpret_cast(rc_param); + return openPMD::switchDatasetType( + implementation::datatype_c_to_cxx(dt), *rc, + implementation::pointer_to_vector( + offset, dimensions, invert), + implementation::pointer_to_vector( + extent, dimensions, invert), + data); +} + +int openPMD_RecordComponent_storeChunkSpan( + // in + openPMD_RecordComponent rc_param, openPMD_Datatype dt, int dimensions, + int const *offset, int const *extent, int invert, + // out + openPMD_DynamicMemoryView *memory_view_param) { + auto rc = reinterpret_cast(rc_param); + auto memory_view = + reinterpret_cast(memory_view_param); + *memory_view = new implementation::DynamicMemoryView(); + (*memory_view)->dtype = dt; + return openPMD::switchDatasetType< + implementation::RecordComponent_storeChunkSpan>( + implementation::datatype_c_to_cxx(dt), *rc, + implementation::pointer_to_vector( + offset, dimensions, invert), + implementation::pointer_to_vector( + extent, dimensions, invert), + (*memory_view)->memory_view); +} + +int openPMD_DynamicMemoryView_resolve( + // in + openPMD_DynamicMemoryView memory_view_param, int deallocate, + // out + void **write_buffer) { + auto memory_view = + reinterpret_cast(memory_view_param); + *write_buffer = + openPMD::switchDatasetType( + implementation::datatype_c_to_cxx(memory_view->dtype), *memory_view); + if (deallocate) { + delete memory_view; + } + return 0; +} + +int openPMD_ParticleSpecies_get_Record( + // in + openPMD_ParticleSpecies ps_param, char const *name, + // out + openPMD_Record *record) { + auto ps = reinterpret_cast(ps_param); + auto &res = ps->operator[](name); + *reinterpret_cast(record) = &res; + return 0; +} + +int openPMD_Record_upcast_to_RecordComponent( + // in + openPMD_Record record, + // out + openPMD_RecordComponent *rc) { + return implementation::do_upcast( + record, rc); +} + +int openPMD_Record_get_Component( + // in + openPMD_Record record_param, char const *name, + // out + openPMD_RecordComponent *rc) { + auto record = reinterpret_cast(record_param); + auto &res = record->operator[](name); + *reinterpret_cast(rc) = &res; + return 0; +} + +char *openPMD_json_merge(char const *into, char const *from, + MPI_Comm maybe_comm) { +#if OPENPMDAPI_VERSION_GE(0, 17, 0) + auto const res = maybe_comm ? openPMD::json::merge(into, from, maybe_comm) + : openPMD::json::merge(into, from); +#else + // openPMD::json::merge in release 0.16.1 has an API bug in which it does not + // resolve filenames, so we do it manually. + auto comm = + maybe_comm ? std::make_optional(maybe_comm) : std::nullopt; + auto const into_resolved = implementation::resolveFileContent(into, comm); + auto const from_resolved = implementation::resolveFileContent(from, comm); + auto const res = std::visit( + [&](auto &&i) { + return std::visit( + [&](auto &&f) { + return openPMD::json::merge( + // poor man's std::forward() + static_cast(i), static_cast(f)); + }, + from_resolved); + }, + into_resolved); +#endif + char *c_res = static_cast(malloc(sizeof(char) * (res.size() + 1))); + std::copy_n(res.c_str(), res.size() + 1, c_res); + return c_res; +} +} + +///////////////// +// auxiliary.h // +///////////////// + +// #ifndef CP2K_OPENPMD_AUXILIARY_H +// #define CP2K_OPENPMD_AUXILIARY_H + +#include + +#ifdef __cplusplus +extern "C" { +#endif // defined(__cplusplus) + +/* Wrapper for MPI_Comm_f2c to ensure it's a compiled symbol. */ +void CP2K_MPI_Comm_f2c( + // in + MPI_Fint comm_f, + // out + MPI_Comm *comm_c); + +size_t CP2K_strlen(char const *); + +void CP2K_free(void *); + +#ifdef __cplusplus +} // extern "C" +#endif // defined(__cplusplus) + +// #endif // !defined(CP2K_OPENPMD_AUXILIARY_H) + +///////////////// +// auxiliary.c // +///////////////// + +#include +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif // defined(__cplusplus) + +void CP2K_MPI_Comm_f2c(MPI_Fint comm_f, MPI_Comm *comm_c) { + *comm_c = MPI_Comm_f2c(comm_f); +} + +size_t CP2K_strlen(char const *s) { return strlen(s); } + +void CP2K_free(void *ptr) { free(ptr); } + +#ifdef __cplusplus +} +#endif // defined(__cplusplus) +#endif // defined(__OPENPMD) diff --git a/src/pw/realspace_grid_openpmd.F b/src/pw/realspace_grid_openpmd.F new file mode 100644 index 0000000000..041e011335 --- /dev/null +++ b/src/pw/realspace_grid_openpmd.F @@ -0,0 +1,526 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Generate Gaussian cube files +! ************************************************************************************************** +MODULE realspace_grid_openpmd + +#ifdef __OPENPMD + USE cp_files, ONLY: close_file, & + open_file + USE cp_log_handling, ONLY: cp_logger_get_default_io_unit + USE cp_output_handling_openpmd, ONLY: cp_openpmd_get_value_unit_nr, & + cp_openpmd_per_call_value_type + USE kinds, ONLY: default_string_length, & + dp + USE message_passing, ONLY: & + file_amode_rdonly, file_offset, mp_comm_type, mp_file_descriptor_type, mp_file_type, & + mp_file_type_free, mp_file_type_hindexed_make_chv, mp_file_type_set_view_chv, & + mpi_character_size +#if defined(__parallel) +#if defined(__MPI_F08) + USE mpi_f08, ONLY: mpi_allreduce, mpi_integer, mpi_bxor, mpi_allgather +#else + USE mpi, ONLY: mpi_allreduce, mpi_integer, mpi_bxor, mpi_allgather +#endif +#endif + + USE openpmd_api, ONLY: & + openpmd_attributable_type, & + openpmd_dynamic_memory_view_type_1d, & + openpmd_dynamic_memory_view_type_3d, & + openpmd_mesh_type, & + openpmd_particle_species_type, & + openpmd_record_component_type, & + openpmd_record_type, openpmd_type_double, openpmd_type_int + USE pw_grid_types, ONLY: PW_MODE_LOCAL + USE pw_types, ONLY: pw_r3d_rs_type + USE util, ONLY: sort_unique + +#else + + USE pw_types, ONLY: pw_r3d_rs_type + USE kinds, ONLY: dp + +#endif + +#include "../base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + PUBLIC ::pw_to_openpmd + +#ifdef __OPENPMD + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'realspace_grid_openpmd' + LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. + LOGICAL, PRIVATE :: parses_linebreaks = .FALSE., & + parse_test = .TRUE. + + TYPE cp_openpmd_write_buffer_1d + REAL(KIND=dp), POINTER :: buffer(:) + END TYPE cp_openpmd_write_buffer_1d +#endif + +CONTAINS + +#ifdef __OPENPMD +! ************************************************************************************************** +!> \brief ... +!> \param particles_z ... +!> \param res_atom_types ... +!> \param res_atom_counts ... +!> \param res_len ... +! ************************************************************************************************** + SUBROUTINE pw_get_atom_types(particles_z, res_atom_types, res_atom_counts, res_len) + INTEGER, DIMENSION(:), INTENT(IN) :: particles_z + INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: res_atom_types, res_atom_counts + INTEGER, INTENT(OUT) :: res_len + + INTEGER :: current_atom_number, i + INTEGER, ALLOCATABLE, DIMENSION(:) :: particles_z_sorted + LOGICAL :: unique + + ALLOCATE (particles_z_sorted(SIZE(particles_z))) + particles_z_sorted(:) = particles_z(:) + CALL sort_unique(particles_z_sorted, unique) + + ALLOCATE (res_atom_types(MIN(118, SIZE(particles_z)))) + ALLOCATE (res_atom_counts(MIN(118, SIZE(particles_z)))) + current_atom_number = -1 + res_len = 0 + DO i = 1, SIZE(particles_z_sorted) + IF (particles_z_sorted(i) /= current_atom_number) THEN + res_len = res_len + 1 + current_atom_number = particles_z_sorted(i) + res_atom_types(res_len) = current_atom_number + res_atom_counts(res_len) = 1 + ELSE + res_atom_counts(res_len) = res_atom_counts(res_len) + 1 + END IF + END DO + + END SUBROUTINE pw_get_atom_types + +! ************************************************************************************************** +!> \brief ... +!> \param particles_z ... +!> \param particles_r ... +!> \param particles_zeff ... +!> \param atom_type ... +!> \param atom_count ... +!> \param openpmd_data ... +!> \param do_write_data ... +! ************************************************************************************************** + SUBROUTINE pw_write_particle_species( & + particles_z, & + particles_r, & + particles_zeff, & + atom_type, & + atom_count, & + openpmd_data, & + do_write_data & + ) + INTEGER, DIMENSION(:), INTENT(IN) :: particles_z + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: particles_r + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: particles_zeff + INTEGER, INTENT(IN) :: atom_type, atom_count + TYPE(cp_openpmd_per_call_value_type) :: openpmd_data + LOGICAL :: do_write_data + + CHARACTER(len=1), DIMENSION(3), PARAMETER :: dims = ["x", "y", "z"] + + CHARACTER(len=3) :: atom_type_as_string + CHARACTER(len=default_string_length) :: species_name + INTEGER :: i, j, k + INTEGER, DIMENSION(1) :: global_extent, global_offset, & + local_extent + TYPE(cp_openpmd_write_buffer_1d) :: charge_write_buffer + TYPE(cp_openpmd_write_buffer_1d), DIMENSION(3) :: write_buffers + TYPE(openpmd_attributable_type) :: attr + TYPE(openpmd_dynamic_memory_view_type_1d) :: unresolved_charge_write_buffer + TYPE(openpmd_dynamic_memory_view_type_1d), & + DIMENSION(3) :: unresolved_write_buffers + TYPE(openpmd_particle_species_type) :: species + TYPE(openpmd_record_component_type) :: charge_component, position_component, & + position_offset_component + TYPE(openpmd_record_type) :: charge, position, position_offset + +! TODO: The charge is probably constant per species? +! If yes, we could use a constant component and save storage space + + global_extent(1) = atom_count + IF (do_write_data) THEN + global_offset(1) = 0 + local_extent(1) = atom_count + ELSE + global_offset(1) = 0 + local_extent(1) = 0 + END IF + + WRITE (atom_type_as_string, '(I3)') atom_type + species_name = TRIM(openpmd_data%name_prefix)//"-"//ADJUSTL(atom_type_as_string) + + species = openpmd_data%iteration%get_particle_species(TRIM(species_name)) + + position_offset = species%get_record("positionOffset") + position = species%get_record("position") + DO k = 1, SIZE(dims) + position_offset_component = position_offset%get_component(dims(k)) + CALL position_offset_component%make_constant_zero(openpmd_type_int, global_extent) + position_component = position%get_component(dims(k)) + CALL position_component%reset_dataset(openpmd_type_double, global_extent) + IF (do_write_data) THEN + unresolved_write_buffers(k) = & + position_component%store_chunk_span_1d_double(global_offset, local_extent) + write_buffers(k)%buffer => unresolved_write_buffers(k)%resolve_double(DEALLOCATE=.FALSE.) + END IF + END DO + + IF (PRESENT(particles_zeff)) THEN + charge = species%get_record("charge") + charge_component = charge%as_record_component() + CALL charge_component%reset_dataset(openpmd_type_double, global_extent) + IF (do_write_data) THEN + unresolved_charge_write_buffer = charge_component%store_chunk_span_1d_double(global_offset, local_extent) + charge_write_buffer%buffer => unresolved_charge_write_buffer%resolve_double(DEALLOCATE=.FALSE.) + END IF + END IF + + IF (do_write_data) THEN + ! Resolve Spans for a second time to allow for internal reallocations in BP4 engine of ADIOS2 + DO k = 1, SIZE(dims) + write_buffers(k)%buffer = unresolved_write_buffers(k)%resolve_double(DEALLOCATE=.TRUE.) + END DO + IF (PRESENT(particles_zeff)) THEN + charge_write_buffer%buffer = unresolved_charge_write_buffer%resolve_double(DEALLOCATE=.TRUE.) + END IF + j = 1 + DO i = 1, SIZE(particles_z) + IF (particles_z(i) == atom_type) THEN + DO k = 1, 3 + write_buffers(k)%buffer(j) = particles_r(k, i) + END DO + IF (PRESENT(particles_zeff)) THEN + charge_write_buffer%buffer(j) = particles_zeff(i) + END IF + j = j + 1 + END IF + END DO + END IF + attr = openpmd_data%iteration%as_attributable() + CALL attr%series_flush("hdf5.independent_stores = true") + END SUBROUTINE pw_write_particle_species + +! ************************************************************************************************** +!> \brief ... +!> \param particles_z ... +!> \param particles_r ... +!> \param particles_zeff ... +!> \param atom_types ... +!> \param atom_counts ... +!> \param num_atom_types ... +!> \param openpmd_data ... +!> \param gid ... +! ************************************************************************************************** + SUBROUTINE pw_write_particles( & + particles_z, & + particles_r, & + particles_zeff, & + atom_types, & + atom_counts, & + num_atom_types, & + openpmd_data, & + gid & + ) + INTEGER, DIMENSION(:), INTENT(IN) :: particles_z + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: particles_r + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: particles_zeff + INTEGER, DIMENSION(:), INTENT(IN) :: atom_types, atom_counts + INTEGER, INTENT(IN), TARGET :: num_atom_types + TYPE(cp_openpmd_per_call_value_type) :: openpmd_data + TYPE(mp_comm_type), OPTIONAL :: gid + + INTEGER :: i, mpi_rank + LOGICAL :: do_write_data + + IF (PRESENT(gid)) THEN + CALL gid%get_rank(mpi_rank) + do_write_data = mpi_rank == 0 + ELSE + do_write_data = .TRUE. + END IF + DO i = 1, num_atom_types + CALL pw_write_particle_species( & + particles_z, & + particles_r, & + particles_zeff, & + atom_types(i), & + atom_counts(i), & + openpmd_data, & + do_write_data & + ) + END DO + END SUBROUTINE pw_write_particles + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param unit_nr ... +!> \param title ... +!> \param particles_r ... +!> \param particles_z ... +!> \param particles_zeff ... +!> \param stride ... +!> \param zero_tails ... +!> \param silent ... +!> \param mpi_io ... +! ************************************************************************************************** + SUBROUTINE pw_to_openpmd( & + pw, & + unit_nr, & + title, & + particles_r, & + particles_z, & + particles_zeff, & + stride, & + zero_tails, & + silent, & + mpi_io & + ) + TYPE(pw_r3d_rs_type), INTENT(IN) :: pw + INTEGER :: unit_nr + CHARACTER(*), INTENT(IN), OPTIONAL :: title + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & + OPTIONAL :: particles_r + INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: particles_z + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: particles_zeff + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: stride + LOGICAL, INTENT(IN), OPTIONAL :: zero_tails, silent, mpi_io + + CHARACTER(len=*), PARAMETER :: routineN = 'pw_to_openpmd' + INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6 + + INTEGER :: count1, count2, count3, handle, i, I1, & + I2, I3, iat, L1, L2, L3, my_rank, & + my_stride(3), np, num_atom_types, & + num_pe, U1, U2, U3 + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_counts, atom_types + INTEGER, DIMENSION(3) :: global_extent, local_extent, offset + LOGICAL :: be_silent, my_zero_tails, parallel_write + REAL(KIND=dp), DIMENSION(3) :: grid_spacing + REAL(KIND=dp), POINTER :: write_buffer(:, :, :) + TYPE(cp_openpmd_per_call_value_type) :: openpmd_data + TYPE(mp_comm_type) :: gid + TYPE(openpmd_attributable_type) :: attr + TYPE(openpmd_dynamic_memory_view_type_3d) :: unresolved_write_buffer + TYPE(openpmd_mesh_type) :: mesh + TYPE(openpmd_record_component_type) :: scalar_mesh + + CALL timeset(routineN, handle) + + my_zero_tails = .FALSE. + be_silent = .FALSE. + parallel_write = .FALSE. + gid = pw%pw_grid%para%group + IF (PRESENT(zero_tails)) my_zero_tails = zero_tails + IF (PRESENT(silent)) be_silent = silent + IF (PRESENT(mpi_io)) parallel_write = mpi_io + my_stride = 1 + IF (PRESENT(stride)) THEN + IF (SIZE(stride) /= 1 .AND. SIZE(stride) /= 3) & + CALL cp_abort(__LOCATION__, "STRIDE keyword can accept only 1 "// & + "(the same for X,Y,Z) or 3 values. Correct your input file.") + IF (SIZE(stride) == 1) THEN + DO i = 1, 3 + my_stride(i) = stride(1) + END DO + ELSE + my_stride = stride(1:3) + END IF + CPASSERT(my_stride(1) > 0) + CPASSERT(my_stride(2) > 0) + CPASSERT(my_stride(3) > 0) + END IF + + openpmd_data = cp_openpmd_get_value_unit_nr(unit_nr) + + CPASSERT(PRESENT(particles_z) .EQV. PRESENT(particles_r)) + np = 0 + IF (PRESENT(particles_z)) THEN + CALL pw_get_atom_types(particles_z, atom_types, atom_counts, num_atom_types) + CPASSERT(SIZE(particles_z) == SIZE(particles_r, dim=2)) + np = SIZE(particles_z) + END IF + + DO i = 1, 3 + ! Notes: + ! 1. This loses information on the rotation of the mesh, the mesh is stored + ! without reference to a global coordinate system + ! 2. This assumes that the coordinate system is not sheared + grid_spacing(i) = SQRT(SUM(pw%pw_grid%dh(:, i)**2))*REAL(my_stride(i), dp) + END DO + + IF (PRESENT(particles_z)) THEN + IF (parallel_write) THEN + CALL pw_write_particles( & + particles_z, & + particles_r, & + particles_zeff, & + atom_types, & + atom_counts, & + num_atom_types, & + openpmd_data, & + gid & + ) + ELSE + CALL pw_write_particles( & + particles_z, & + particles_r, & + particles_zeff, & + atom_types, & + atom_counts, & + num_atom_types, & + openpmd_data & + ) + END IF + END IF + + DO iat = 1, 3 + global_extent(iat) = (pw%pw_grid%npts(iat) + my_stride(iat) - 1)/my_stride(iat) + ! '+ 1' because upper end is inclusive, '- 1' for upper gaussian bracket + offset(iat) = ((pw%pw_grid%bounds_local(1, iat) - pw%pw_grid%bounds(1, iat) + my_stride(iat) - 1)/my_stride(iat)) + ! refer local_extent to the global offset first in order to have consistent rounding + local_extent(iat) = ((pw%pw_grid%bounds_local(2, iat) + 1 - pw%pw_grid%bounds(1, iat) + my_stride(iat) - 1)/my_stride(iat)) + END DO + local_extent = local_extent - offset + + mesh = openpmd_data%iteration%get_mesh(TRIM(openpmd_data%name_prefix)) + CALL mesh%set_axis_labels(["x", "y", "z"]) + CALL mesh%set_position([0.5_dp, 0.5_dp, 0.5_dp]) + CALL mesh%set_grid_global_offset([0._dp, 0._dp, 0._dp]) + CALL mesh%set_grid_spacing(grid_spacing) + scalar_mesh = mesh%as_record_component() + CALL scalar_mesh%reset_dataset(openpmd_type_double, global_extent) + + ! shortcut + ! need to adjust L1/U1 for uneven distributions across MPI ranks + ! (when working with a stride, we might have to skip the first n values) + ! so keep this consistent with the offset and local_extent computed above + ! L1 = pw%pw_grid%bounds_local(1, 1) + L1 = offset(1)*my_stride(1) + L2 = pw%pw_grid%bounds_local(1, 2) + L3 = pw%pw_grid%bounds_local(1, 3) + ! U1 = pw%pw_grid%bounds_local(2, 1) + ! offset + local_extent is the start index for the next rank already + ! since the indexes are inclusive, subtract 1 from the boundary index + U1 = (offset(1) + local_extent(1) - 1)*my_stride(1) + U2 = pw%pw_grid%bounds_local(2, 2) + U3 = pw%pw_grid%bounds_local(2, 3) + + my_rank = pw%pw_grid%para%group%mepos + num_pe = pw%pw_grid%para%group%num_pe + + IF (ALL(my_stride == 1)) THEN + CALL scalar_mesh%store_chunk(pw%array(L1:U1, L2:U2, L3:U3), offset) + ! Are there some conditions under which we can skip this flush? + attr = openpmd_data%iteration%as_attributable() + CALL attr%series_flush("hdf5.independent_stores = false") + ELSE + count3 = 0 + DO I3 = L3, U3, my_stride(3) + ! maybe add an overload to provide `buf` here for HDF5, might have better performance + ! for intermittent flushing + ! or just call the buffer in the outer function if memory is no problem... + unresolved_write_buffer = scalar_mesh%store_chunk_span_3d_double( & + [offset(1), offset(2), offset(3) + count3], & + [local_extent(1), local_extent(2), 1]) + write_buffer => unresolved_write_buffer%resolve_double(DEALLOCATE=.TRUE.) + + ! Sanity checks: ensure buffer is associated and matches expected shape + CPASSERT(ASSOCIATED(write_buffer)) + CPASSERT(SIZE(write_buffer, 1) == local_extent(1)) + CPASSERT(SIZE(write_buffer, 2) == local_extent(2)) + CPASSERT(SIZE(write_buffer, 3) == 1) + + count2 = 0 + DO I2 = L2, U2, my_stride(2) + ! This loop deals with ray (:, count2, count3) of the local subspace + ! The write buffer itself has been allocated for slice (:, :, count3) + count1 = 0 + DO I1 = L1, U1, my_stride(1) + write_buffer(count1 + 1, count2 + 1, 1) = pw%array(I1, I2, I3) + ! Debug: print the target indices in write_buffer to the command line + ! WRITE(*,*) 'write_buffer index:', count1 + 1, ',', count2 + 1, ',', 1 + count1 = count1 + 1 + END DO + count2 = count2 + 1 + END DO + count3 = count3 + 1 + END DO + END IF + + CALL timestop(handle) + + END SUBROUTINE pw_to_openpmd + +#else + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param unit_nr ... +!> \param title ... +!> \param particles_r ... +!> \param particles_z ... +!> \param particles_zeff ... +!> \param stride ... +!> \param zero_tails ... +!> \param silent ... +!> \param mpi_io ... +! ************************************************************************************************** + SUBROUTINE pw_to_openpmd( & + pw, & + unit_nr, & + title, & + particles_r, & + particles_z, & + particles_zeff, & + stride, & + zero_tails, & + silent, & + mpi_io & + ) + TYPE(pw_r3d_rs_type), INTENT(IN) :: pw + INTEGER :: unit_nr + CHARACTER(*), INTENT(IN), OPTIONAL :: title + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & + OPTIONAL :: particles_r + INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: particles_z + REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: particles_zeff + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: stride + LOGICAL, INTENT(IN), OPTIONAL :: zero_tails, silent, mpi_io + + MARK_USED(pw) + MARK_USED(unit_nr) + MARK_USED(title) + MARK_USED(particles_r) + MARK_USED(particles_z) + MARK_USED(particles_zeff) + MARK_USED(stride) + MARK_USED(zero_tails) + MARK_USED(silent) + MARK_USED(mpi_io) + CPABORT("CP2K compiled without the openPMD-api") + + END SUBROUTINE pw_to_openpmd + +#endif + +END MODULE realspace_grid_openpmd diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 6f8b190b05..4fb124f3cc 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -51,7 +51,11 @@ MODULE qs_scf_post_gpw cp_print_key_finished_output,& cp_print_key_should_output,& cp_print_key_unit_nr + USE cp_output_handling_openpmd, ONLY: cp_openpmd_close_iterations,& + cp_openpmd_print_key_finished_output,& + cp_openpmd_print_key_unit_nr USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube + USE cp_realspace_grid_openpmd, ONLY: cp_pw_to_openpmd USE dct, ONLY: pw_shrink USE ed_analysis, ONLY: edmf_analysis USE eeq_method, ONLY: eeq_print @@ -215,8 +219,306 @@ MODULE qs_scf_post_gpw PUBLIC :: make_lumo_gpw + CHARACTER(len=*), PARAMETER :: & + str_mo_cubes = "PRINT%MO_CUBES", & + str_mo_openpmd = "PRINT%MO_OPENPMD", & + str_elf_cubes = "PRINT%ELF_CUBE", & + str_elf_openpmd = "PRINT%ELF_OPENPMD", & + str_e_density_cubes = "PRINT%E_DENSITY_CUBE", & + str_e_density_openpmd = "PRINT%E_DENSITY_OPENPMD" + + INTEGER, PARAMETER :: grid_output_cubes = 1, grid_output_openpmd = 2 + + ! Generic information on whether a certain output section has been activated + ! or not, and on whether it has been activated in the Cube or openPMD variant. + ! Create with function cube_or_openpmd(), see there for further details. + TYPE cp_section_key + CHARACTER(len=default_string_length) :: relative_section_key ! e.g. PRINT%MO_CUBES + CHARACTER(len=default_string_length) :: absolute_section_key ! e.g. DFT%PRINT%MO_CUBES + CHARACTER(len=7) :: format_name ! 'openPMD' or 'Cube', for logging + INTEGER :: grid_output ! either 1 for grid_output_cubes or 2 for grid_output_openpmd + LOGICAL :: do_output + CONTAINS + ! Open a file as either Cube or openPMD + PROCEDURE, PUBLIC :: print_key_unit_nr => cp_forward_print_key_unit_nr + ! Write either to the Cube or openPMD file + PROCEDURE, PUBLIC :: write_pw => cp_forward_write_pw + ! Close either the Cube or openPMD file + PROCEDURE, PUBLIC :: print_key_finished_output => cp_forward_print_key_finished_output + ! Helpers + PROCEDURE, PUBLIC :: do_openpmd => cp_section_key_do_openpmd + PROCEDURE, PUBLIC :: do_cubes => cp_section_key_do_cubes + PROCEDURE, PUBLIC :: concat_to_relative => cp_section_key_concat_to_relative + PROCEDURE, PUBLIC :: concat_to_absolute => cp_section_key_concat_to_absolute + END TYPE cp_section_key + CONTAINS +! ************************************************************************************************** +!> \brief Append `extend_by` to the absolute path of the base section. +!> \param self ... +!> \param extend_by ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_section_key_concat_to_absolute(self, extend_by) RESULT(res) + CLASS(cp_section_key), INTENT(IN) :: self + CHARACTER(*), INTENT(IN) :: extend_by + CHARACTER(len=default_string_length) :: res + + IF (LEN(TRIM(extend_by)) > 0 .AND. extend_by(1:1) == "%") THEN + res = TRIM(self%absolute_section_key)//TRIM(extend_by) + ELSE + res = TRIM(self%absolute_section_key)//"%"//TRIM(extend_by) + END IF + END FUNCTION cp_section_key_concat_to_absolute + +! ************************************************************************************************** +!> \brief Append `extend_by` to the relative path (e.g. without DFT%) of the base section. +!> \param self ... +!> \param extend_by ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_section_key_concat_to_relative(self, extend_by) RESULT(res) + CLASS(cp_section_key), INTENT(IN) :: self + CHARACTER(*), INTENT(IN) :: extend_by + CHARACTER(len=default_string_length) :: res + + IF (LEN(TRIM(extend_by)) > 0 .AND. extend_by(1:1) == "%") THEN + res = TRIM(self%relative_section_key)//TRIM(extend_by) + ELSE + res = TRIM(self%relative_section_key)//"%"//TRIM(extend_by) + END IF + END FUNCTION cp_section_key_concat_to_relative + +! ************************************************************************************************** +!> \brief Is Cube output active for the current base section? +!> \param self ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_section_key_do_cubes(self) RESULT(res) + CLASS(cp_section_key) :: self + LOGICAL :: res + + res = self%do_output .AND. self%grid_output == grid_output_cubes + END FUNCTION cp_section_key_do_cubes + +! ************************************************************************************************** +!> \brief Is openPMD output active for the current base section? +!> \param self ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_section_key_do_openpmd(self) RESULT(res) + CLASS(cp_section_key) :: self + LOGICAL :: res + + res = self%do_output .AND. self%grid_output == grid_output_openpmd + END FUNCTION cp_section_key_do_openpmd + +! ************************************************************************************************** +!> \brief Forwards to either `cp_print_key_unit_nr` or `cp_openpmd_print_key_unit_nr`, +!> depending on the configuration of the current base section. +!> Opens either a Cube or openPMD output file +!> \param self ... +!> \param logger ... +!> \param basis_section ... +!> \param print_key_path ... +!> \param extension ... +!> \param middle_name ... +!> \param local ... +!> \param log_filename ... +!> \param ignore_should_output ... +!> \param file_form ... +!> \param file_position ... +!> \param file_action ... +!> \param file_status ... +!> \param do_backup ... +!> \param on_file ... +!> \param is_new_file ... +!> \param mpi_io ... +!> \param fout ... +!> \param openpmd_basename ... +!> \return ... +! ************************************************************************************************** + FUNCTION cp_forward_print_key_unit_nr(self, logger, basis_section, print_key_path, extension, & + middle_name, local, log_filename, ignore_should_output, file_form, file_position, & + file_action, file_status, do_backup, on_file, is_new_file, mpi_io, & + fout, openpmd_basename) RESULT(res) + CLASS(cp_section_key), INTENT(IN) :: self + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), INTENT(IN) :: basis_section + CHARACTER(len=*), INTENT(IN), OPTIONAL :: print_key_path + CHARACTER(len=*), INTENT(IN) :: extension + CHARACTER(len=*), INTENT(IN), OPTIONAL :: middle_name + LOGICAL, INTENT(IN), OPTIONAL :: local, log_filename, ignore_should_output + CHARACTER(len=*), INTENT(IN), OPTIONAL :: file_form, file_position, file_action, & + file_status + LOGICAL, INTENT(IN), OPTIONAL :: do_backup, on_file + LOGICAL, INTENT(OUT), OPTIONAL :: is_new_file + LOGICAL, INTENT(INOUT), OPTIONAL :: mpi_io + CHARACTER(len=default_path_length), INTENT(OUT), & + OPTIONAL :: fout + CHARACTER(len=*), INTENT(IN), OPTIONAL :: openpmd_basename + INTEGER :: res + + IF (self%grid_output == grid_output_cubes) THEN + res = cp_print_key_unit_nr( & + logger, basis_section, print_key_path, extension=extension, & + middle_name=middle_name, local=local, log_filename=log_filename, & + ignore_should_output=ignore_should_output, file_form=file_form, & + file_position=file_position, file_action=file_action, & + file_status=file_status, do_backup=do_backup, on_file=on_file, & + is_new_file=is_new_file, mpi_io=mpi_io, fout=fout) + ELSE + res = cp_openpmd_print_key_unit_nr(logger, basis_section, print_key_path, & + middle_name=middle_name, ignore_should_output=ignore_should_output, & + mpi_io=mpi_io, fout=fout, openpmd_basename=openpmd_basename) + END IF + END FUNCTION cp_forward_print_key_unit_nr + +! ************************************************************************************************** +!> \brief Forwards to either `cp_pw_to_cube` or `cp_pw_to_openpmd`, +!> depending on the configuration of the current base section. +!> Writes data to either a Cube or an openPMD file. +!> \param self ... +!> \param pw ... +!> \param unit_nr ... +!> \param title ... +!> \param particles ... +!> \param zeff ... +!> \param stride ... +!> \param max_file_size_mb ... +!> \param zero_tails ... +!> \param silent ... +!> \param mpi_io ... +! ************************************************************************************************** + SUBROUTINE cp_forward_write_pw( & + self, & + pw, & + unit_nr, & + title, & + particles, & + zeff, & + stride, & + max_file_size_mb, & + zero_tails, & + silent, & + mpi_io & + ) + CLASS(cp_section_key), INTENT(IN) :: self + TYPE(pw_r3d_rs_type), INTENT(IN) :: pw + INTEGER, INTENT(IN) :: unit_nr + CHARACTER(*), INTENT(IN), OPTIONAL :: title + TYPE(particle_list_type), POINTER :: particles + INTEGER, DIMENSION(:), OPTIONAL, POINTER :: stride + REAL(KIND=dp), INTENT(IN), OPTIONAL :: max_file_size_mb + LOGICAL, INTENT(IN), OPTIONAL :: zero_tails, silent, mpi_io + REAL(KIND=dp), DIMENSION(:), OPTIONAL :: zeff + + IF (self%grid_output == grid_output_cubes) THEN + CALL cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io) + ELSE + CALL cp_pw_to_openpmd(pw, unit_nr, title, particles, zeff, stride, zero_tails, silent, mpi_io) + END IF + END SUBROUTINE cp_forward_write_pw + +! ************************************************************************************************** +!> \brief Forwards to either `cp_print_key_finished_output` or `cp_openpmd_print_key_finished_output`, +!> depending on the configuration of the current base section. +!> Closes either a Cube file or a reference to a section within an openPMD file. +!> \param self ... +!> \param unit_nr ... +!> \param logger ... +!> \param basis_section ... +!> \param print_key_path ... +!> \param local ... +!> \param ignore_should_output ... +!> \param on_file ... +!> \param mpi_io ... +! ************************************************************************************************** + SUBROUTINE cp_forward_print_key_finished_output(self, unit_nr, logger, basis_section, & + print_key_path, local, ignore_should_output, on_file, & + mpi_io) + CLASS(cp_section_key), INTENT(IN) :: self + INTEGER, INTENT(INOUT) :: unit_nr + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), INTENT(IN) :: basis_section + CHARACTER(len=*), INTENT(IN), OPTIONAL :: print_key_path + LOGICAL, INTENT(IN), OPTIONAL :: local, ignore_should_output, on_file, & + mpi_io + + IF (self%grid_output == grid_output_cubes) THEN + CALL cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io) + ELSE + CALL cp_openpmd_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, mpi_io) + END IF + END SUBROUTINE cp_forward_print_key_finished_output + + ! +! ************************************************************************************************** +!> \brief Decides if a particular output routine will write to openPMD, to Cube or to none. +!> Writing to both is not supported. +!> The distinction between Cube and openPMD output works such that the output configuration +!> sections exist as duplicates: E.g. for DFT%PRINT%MO_CUBES, +!> there additionally exists DFT%PRINT%MO_OPENPMD. +!> The internal base configuration for such sections is identical; additionally there +!> exist format-specific options such as APPEND for Cube or OPENPMD_CFG_FILE for openPMD. +!> The routines in this file alternate between using relative section paths without the +!> %DFT prefix (e.g. PRINT%MO_CUBES) or absolute section paths with the %DF% prefix +!> (e.g. DFT%PRINT%MO_CUBES). Call this routine with the relative paths. +!> \param input ... +!> \param str_cubes ... +!> \param str_openpmd ... +!> \param logger ... +!> \return ... +! ************************************************************************************************** + FUNCTION cube_or_openpmd(input, str_cubes, str_openpmd, logger) RESULT(res) + TYPE(section_vals_type), POINTER :: input + CHARACTER(len=*), INTENT(IN) :: str_cubes, str_openpmd + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_section_key) :: res + + LOGICAL :: do_cubes, do_openpmd + + do_cubes = BTEST(cp_print_key_should_output( & + logger%iter_info, input, & + "DFT%"//TRIM(ADJUSTL(str_cubes))), cp_p_file) + do_openpmd = BTEST(cp_print_key_should_output( & + logger%iter_info, input, & + "DFT%"//TRIM(ADJUSTL(str_openpmd))), cp_p_file) + ! Having Cube and openPMD output both active should be theoretically possible. + ! It would require some extra handling for the unit_nr return values. + ! (e.g. returning the Cube unit_nr and internally storing the associated openPMD unit_nr). + CPASSERT(.NOT. (do_cubes .AND. do_openpmd)) + res%do_output = do_cubes .OR. do_openpmd + IF (do_openpmd) THEN + res%grid_output = grid_output_openpmd + res%relative_section_key = TRIM(ADJUSTL(str_openpmd)) + res%format_name = "openPMD" + ELSE + res%grid_output = grid_output_cubes + res%relative_section_key = TRIM(ADJUSTL(str_cubes)) + res%format_name = "Cube" + END IF + res%absolute_section_key = "DFT%"//TRIM(ADJUSTL(res%relative_section_key)) + END FUNCTION cube_or_openpmd + +! ************************************************************************************************** +!> \brief This section key is named WRITE_CUBE for Cube which does not make much sense +!> for openPMD, so this key name has to be distinguished. +!> \param grid_output ... +!> \return ... +! ************************************************************************************************** + FUNCTION section_key_do_write(grid_output) RESULT(res) + INTEGER, INTENT(IN) :: grid_output + CHARACTER(len=32) :: res + + IF (grid_output == grid_output_cubes) THEN + res = "%WRITE_CUBE" + ELSE IF (grid_output == grid_output_openpmd) THEN + res = "%WRITE_OPENPMD" + END IF + END FUNCTION section_key_do_write + ! ************************************************************************************************** !> \brief collects possible post - scf calculations and prints info / computes properties. !> \param qs_env the qs_env in which the qs_env lives @@ -241,14 +543,16 @@ CONTAINS CHARACTER(6), OPTIONAL :: wf_type LOGICAL, OPTIONAL :: do_mp2 - CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw' + CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw', & + warning_cube_kpoint = "Print MO cubes not implemented for k-point calculations", & + warning_openpmd_kpoint = "Writing to openPMD not implemented for k-point calculations" INTEGER :: handle, homo, ispin, min_lumos, n_rep, & nchk_nmoloc, nhomo, nlumo, nlumo_stm, & nlumos, nmo, nspins, output_unit, & unit_nr INTEGER, DIMENSION(:, :, :), POINTER :: marked_states - LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mixed, do_mo_cubes, do_stm, & + LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mixed, do_stm, & do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_do_mp2, & my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo, p_loc_mixed REAL(dp) :: e_kin @@ -266,6 +570,7 @@ CONTAINS unoccupied_orbs, unoccupied_orbs_stm TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_section_key) :: mo_section TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_p_mp2, matrix_s, & mo_derivs TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: kinetic_m, rho_ao @@ -435,20 +740,20 @@ CONTAINS nlumo_stm = 0 IF (do_stm) nlumo_stm = section_get_ival(stm_section, "NLUMO") - ! check for CUBES (MOs and WANNIERS) - do_mo_cubes = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") & - , cp_p_file) + ! check for CUBES or openPMD (MOs and WANNIERS) + mo_section = cube_or_openpmd(input, str_mo_cubes, str_mo_openpmd, logger) + IF (loc_print_explicit) THEN do_wannier_cubes = BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, & "WANNIER_CUBES"), cp_p_file) ELSE do_wannier_cubes = .FALSE. END IF - nlumo = section_get_ival(dft_section, "PRINT%MO_CUBES%NLUMO") - nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO") + nlumo = section_get_ival(dft_section, mo_section%concat_to_relative("%NLUMO")) + nhomo = section_get_ival(dft_section, mo_section%concat_to_relative("%NHOMO")) ! Setup the grids needed to compute a wavefunction given a vector.. - IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN + IF (((mo_section%do_output .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & pw_pools=pw_pools) CALL auxbas_pw_pool%create_pw(wf_r) @@ -462,18 +767,19 @@ CONTAINS nspins = dft_control%nspins END IF !Some info about ROKS - IF (dft_control%restricted .AND. (do_mo_cubes .OR. p_loc_homo)) THEN + IF (dft_control%restricted .AND. (mo_section%do_output .OR. p_loc_homo)) THEN CALL cp_abort(__LOCATION__, "Unclear how we define MOs / localization in the restricted case ... ") ! It is possible to obtain Wannier centers for ROKS without rotations for SINGLE OCCUPIED ORBITALS END IF ! Makes the MOs eigenstates, computes eigenvalues, write cubes IF (do_kpoints) THEN - CPWARN_IF(do_mo_cubes, "Print MO cubes not implemented for k-point calculations") + CPWARN_IF(mo_section%do_cubes(), warning_cube_kpoint) + CPWARN_IF(mo_section%do_openpmd(), warning_openpmd_kpoint) ELSE CALL get_qs_env(qs_env, & mos=mos, & matrix_ks=ks_rmpv) - IF ((do_mo_cubes .AND. nhomo /= 0) .OR. do_stm) THEN + IF ((mo_section%do_output .AND. nhomo /= 0) .OR. do_stm) THEN CALL get_qs_env(qs_env, mo_derivs=mo_derivs) IF (dft_control%do_admm) THEN CALL get_qs_env(qs_env, admm_env=admm_env) @@ -494,13 +800,13 @@ CONTAINS END DO has_homo = .TRUE. END IF - IF (do_mo_cubes .AND. nhomo /= 0) THEN + IF (mo_section%do_output .AND. nhomo /= 0) THEN DO ispin = 1, nspins ! Prints the cube files of OCCUPIED ORBITALS CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, & eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo) CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, & - mo_coeff, wf_g, wf_r, particles, homo, ispin) + mo_coeff, wf_g, wf_r, particles, homo, ispin, mo_section) END DO END IF END IF @@ -536,7 +842,7 @@ CONTAINS CALL qs_loc_env_create(qs_loc_env_homo) CALL qs_loc_control_init(qs_loc_env_homo, localize_section, do_homo=do_homo) CALL qs_loc_init(qs_env, qs_loc_env_homo, localize_section, homo_localized, do_homo, & - do_mo_cubes, mo_loc_history=mo_loc_history) + mo_section%do_output, mo_loc_history=mo_loc_history) CALL get_localization_info(qs_env, qs_loc_env_homo, localize_section, homo_localized, & wf_r, wf_g, particles, occupied_orbs, occupied_evals, marked_states) @@ -561,12 +867,12 @@ CONTAINS ! Gets the lumos, and eigenvalues for the lumos, and localize them if requested IF (do_kpoints) THEN - IF (do_mo_cubes .OR. p_loc_lumo) THEN + IF (mo_section%do_output .OR. p_loc_lumo) THEN ! nothing at the moment, not implemented CPWARN("Localization and MO related output not implemented for k-point calculations!") END IF ELSE - compute_lumos = do_mo_cubes .AND. nlumo /= 0 + compute_lumos = mo_section%do_output .AND. nlumo /= 0 compute_lumos = compute_lumos .OR. p_loc_lumo DO ispin = 1, dft_control%nspins @@ -574,9 +880,9 @@ CONTAINS compute_lumos = compute_lumos .AND. homo == nmo END DO - IF (do_mo_cubes .AND. .NOT. compute_lumos) THEN + IF (mo_section%do_output .AND. .NOT. compute_lumos) THEN - nlumo = section_get_ival(dft_section, "PRINT%MO_CUBES%NLUMO") + nlumo = section_get_ival(dft_section, mo_section%concat_to_relative("%NLUMO")) DO ispin = 1, dft_control%nspins CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo, eigenvalues=mo_eigenvalues) @@ -594,7 +900,7 @@ CONTAINS ! Prints the cube files of UNOCCUPIED ORBITALS CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff) CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & - mo_coeff, wf_g, wf_r, particles, nlumo, homo, ispin, lumo=homo + 1) + mo_coeff, wf_g, wf_r, particles, nlumo, homo, ispin, lumo=homo + 1, mo_section=mo_section) END IF END DO @@ -624,7 +930,7 @@ CONTAINS IF (p_loc_lumo .AND. nlumo /= -1) nlumos = MIN(nlumo, nlumos) ! Prints the cube files of UNOCCUPIED ORBITALS CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & - unoccupied_orbs(ispin), wf_g, wf_r, particles, nlumos, homo, ispin) + unoccupied_orbs(ispin), wf_g, wf_r, particles, nlumos, homo, ispin, mo_section=mo_section) END IF END DO @@ -634,7 +940,7 @@ CONTAINS CALL cp_fm_create(lumo_localized(ispin), unoccupied_orbs(ispin)%matrix_struct) CALL cp_fm_to_fm(unoccupied_orbs(ispin), lumo_localized(ispin)) END DO - CALL qs_loc_init(qs_env, qs_loc_env_lumo, localize_section, lumo_localized, do_homo, do_mo_cubes, & + CALL qs_loc_init(qs_env, qs_loc_env_lumo, localize_section, lumo_localized, do_homo, mo_section%do_output, & evals=unoccupied_evals) CALL qs_loc_env_init(qs_loc_env_lumo, qs_loc_env_lumo%localized_wfn_control, qs_env, & loc_coeff=unoccupied_orbs) @@ -687,7 +993,7 @@ CONTAINS CALL qs_loc_env_create(qs_loc_env_mixed) CALL qs_loc_control_init(qs_loc_env_mixed, localize_section, do_homo=do_homo, do_mixed=do_mixed) CALL qs_loc_init(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, do_homo, & - do_mo_cubes, mo_loc_history=mo_loc_history, tot_zeff_corr=total_zeff_corr, & + mo_section%do_output, mo_loc_history=mo_loc_history, tot_zeff_corr=total_zeff_corr, & do_mixed=do_mixed) DO ispin = 1, dft_control%nspins @@ -718,7 +1024,7 @@ CONTAINS END IF ! Deallocate grids needed to compute wavefunctions - IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN + IF (((mo_section%do_output .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN CALL auxbas_pw_pool%give_back_pw(wf_r) CALL auxbas_pw_pool%give_back_pw(wf_g) END IF @@ -822,6 +1128,8 @@ CONTAINS CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) END IF + CALL cp_openpmd_close_iterations() + CALL timestop(handle) END SUBROUTINE scf_post_calculation_gpw @@ -1001,9 +1309,10 @@ CONTAINS !> \param particles ... !> \param homo ... !> \param ispin ... +!> \param mo_section ... ! ************************************************************************************************** SUBROUTINE qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, & - mo_coeff, wf_g, wf_r, particles, homo, ispin) + mo_coeff, wf_g, wf_r, particles, homo, ispin, mo_section) TYPE(section_vals_type), POINTER :: input, dft_section TYPE(dft_control_type), POINTER :: dft_control TYPE(cp_logger_type), POINTER :: logger @@ -1013,6 +1322,7 @@ CONTAINS TYPE(pw_r3d_rs_type), INTENT(INOUT) :: wf_r TYPE(particle_list_type), POINTER :: particles INTEGER, INTENT(IN) :: homo, ispin + TYPE(cp_section_key) :: mo_section CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_occ_cubes' @@ -1029,22 +1339,31 @@ CONTAINS CALL timeset(routineN, handle) +#ifndef __OPENPMD + ! Error should usually be caught earlier as PRINT%MO_OPENPMD is not added to the input section + ! if openPMD is not activated + CPASSERT(mo_section%grid_output /= grid_output_openpmd) +#endif + NULLIFY (list_index) - IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") & - , cp_p_file) .AND. section_get_lval(dft_section, "PRINT%MO_CUBES%WRITE_CUBE")) THEN - nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO") - append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND") + IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, mo_section%relative_section_key) & + , cp_p_file) .AND. section_get_lval(dft_section, mo_section%concat_to_relative(section_key_do_write(mo_section%grid_output)))) THEN + nhomo = section_get_ival(dft_section, mo_section%concat_to_relative("%NHOMO")) + ! For openPMD, refer to access modes instead of APPEND key + IF (mo_section%grid_output == grid_output_cubes) THEN + append_cube = section_get_lval(dft_section, mo_section%concat_to_relative("%APPEND")) + END IF my_pos_cube = "REWIND" IF (append_cube) THEN my_pos_cube = "APPEND" END IF - CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", n_rep_val=n_rep) + CALL section_vals_val_get(dft_section, mo_section%concat_to_relative("%HOMO_LIST"), n_rep_val=n_rep) IF (n_rep > 0) THEN ! write the cubes of the list nlist = 0 DO ir = 1, n_rep NULLIFY (list) - CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", i_rep_val=ir, & + CALL section_vals_val_get(dft_section, mo_section%concat_to_relative("%HOMO_LIST"), i_rep_val=ir, & i_vals=list) IF (ASSOCIATED(list)) THEN CALL reallocate(list_index, 1, nlist + SIZE(list)) @@ -1075,15 +1394,16 @@ CONTAINS cell, dft_control, particle_set, pw_env) WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", & - middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & - mpi_io=mpi_io) + + unit_nr = mo_section%print_key_unit_nr(logger, input, mo_section%absolute_section_key, extension=".cube", & + middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & + mpi_io=mpi_io, openpmd_basename="dft-mo") WRITE (title, *) "WAVEFUNCTION ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo - CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, & - stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), & - max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), & - mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io) + CALL mo_section%write_pw(wf_r, unit_nr, title, particles=particles, & + stride=section_get_ivals(dft_section, mo_section%concat_to_relative("%STRIDE")), & + max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), & + mpi_io=mpi_io) + CALL mo_section%print_key_finished_output(unit_nr, logger, input, mo_section%absolute_section_key, mpi_io=mpi_io) END DO IF (ASSOCIATED(list_index)) DEALLOCATE (list_index) END IF @@ -1107,9 +1427,10 @@ CONTAINS !> \param homo ... !> \param ispin ... !> \param lumo ... +!> \param mo_section ... ! ************************************************************************************************** SUBROUTINE qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & - unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, lumo) + unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, lumo, mo_section) TYPE(section_vals_type), POINTER :: input, dft_section TYPE(dft_control_type), POINTER :: dft_control @@ -1121,6 +1442,7 @@ CONTAINS TYPE(particle_list_type), POINTER :: particles INTEGER, INTENT(IN) :: nlumos, homo, ispin INTEGER, INTENT(IN), OPTIONAL :: lumo + TYPE(cp_section_key) :: mo_section CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_unocc_cubes' @@ -1136,10 +1458,19 @@ CONTAINS CALL timeset(routineN, handle) - IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES"), cp_p_file) & - .AND. section_get_lval(dft_section, "PRINT%MO_CUBES%WRITE_CUBE")) THEN +#ifndef __OPENPMD + ! Error should usually be caught earlier as PRINT%MO_OPENPMD is not added to the input section + ! if openPMD is not activated + CPASSERT(mo_section%grid_output /= grid_output_openpmd) +#endif + + IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, mo_section%relative_section_key), cp_p_file) & + .AND. section_get_lval(dft_section, mo_section%concat_to_relative(section_key_do_write(mo_section%grid_output)))) THEN NULLIFY (qs_kind_set, particle_set, pw_env, cell) - append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND") + ! For openPMD, refer to access modes instead of APPEND key + IF (mo_section%grid_output == grid_output_cubes) THEN + append_cube = section_get_lval(dft_section, mo_section%concat_to_relative("%APPEND")) + END IF my_pos_cube = "REWIND" IF (append_cube) THEN my_pos_cube = "APPEND" @@ -1164,15 +1495,16 @@ CONTAINS WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", index_mo, "_", ispin mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", & - middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & - mpi_io=mpi_io) + unit_nr = mo_section%print_key_unit_nr(logger, input, mo_section%absolute_section_key, extension=".cube", & + middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & + mpi_io=mpi_io, openpmd_basename="dft-mo") WRITE (title, *) "WAVEFUNCTION ", index_mo, " spin ", ispin, " i.e. LUMO + ", ifirst + ivector - 2 - CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, & - stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), & - max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), & - mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io) + CALL mo_section%write_pw(wf_r, unit_nr, title, particles=particles, & + stride=section_get_ivals(dft_section, mo_section%concat_to_relative("%STRIDE")), & + max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), & + mpi_io=mpi_io) + CALL mo_section%print_key_finished_output(unit_nr, logger, input, mo_section%absolute_section_key, mpi_io=mpi_io) + END DO END IF @@ -1448,6 +1780,7 @@ CONTAINS INTEGER :: handle, ispin, output_unit, unit_nr LOGICAL :: append_cube, gapw, mpi_io REAL(dp) :: rho_cutoff + TYPE(cp_section_key) :: elf_section_key TYPE(dft_control_type), POINTER :: dft_control TYPE(particle_list_type), POINTER :: particles TYPE(pw_env_type), POINTER :: pw_env @@ -1460,9 +1793,10 @@ CONTAINS CALL timeset(routineN, handle) output_unit = cp_logger_get_default_io_unit(logger) - elf_section => section_vals_get_subs_vals(input, "DFT%PRINT%ELF_CUBE") - IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & - "DFT%PRINT%ELF_CUBE"), cp_p_file)) THEN + elf_section_key = cube_or_openpmd(input, str_elf_cubes, str_elf_openpmd, logger) + + elf_section => section_vals_get_subs_vals(input, elf_section_key%absolute_section_key) + IF (elf_section_key%do_output) THEN NULLIFY (dft_control, pw_env, auxbas_pw_pool, pw_pools, particles, subsys) CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, subsys=subsys) @@ -1487,7 +1821,11 @@ CONTAINS CALL qs_elf_calc(qs_env, elf_r, rho_cutoff) ! write ELF into cube file - append_cube = section_get_lval(elf_section, "APPEND") + + ! For openPMD, refer to access modes instead of APPEND key + IF (elf_section_key%grid_output == grid_output_cubes) THEN + append_cube = section_get_lval(elf_section, "APPEND") + END IF my_pos_cube = "REWIND" IF (append_cube) THEN my_pos_cube = "APPEND" @@ -1497,9 +1835,10 @@ CONTAINS WRITE (filename, '(a5,I1.1)') "ELF_S", ispin WRITE (title, *) "ELF spin ", ispin mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%ELF_CUBE", extension=".cube", & - middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & - mpi_io=mpi_io, fout=mpi_filename) + unit_nr = elf_section_key%print_key_unit_nr( & + logger, input, elf_section_key%absolute_section_key, extension=".cube", & + middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & + mpi_io=mpi_io, fout=mpi_filename, openpmd_basename="dft-elf") IF (output_unit > 0) THEN IF (.NOT. mpi_io) THEN INQUIRE (UNIT=unit_nr, NAME=filename) @@ -1507,13 +1846,18 @@ CONTAINS filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "ELF is written in cube file format to the file:", & + "ELF is written in "//elf_section_key%format_name//" file format to the file:", & TRIM(filename) END IF - CALL cp_pw_to_cube(elf_r(ispin), unit_nr, title, particles=particles, & - stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%ELF_CUBE", mpi_io=mpi_io) + CALL elf_section_key%write_pw(elf_r(ispin), unit_nr, title, particles=particles, & + stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io) + CALL elf_section_key%print_key_finished_output( & + unit_nr, & + logger, & + input, & + elf_section_key%absolute_section_key, & + mpi_io=mpi_io) CALL auxbas_pw_pool%give_back_pw(elf_r(ispin)) END DO @@ -1905,6 +2249,7 @@ CONTAINS TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_section_key) :: e_density_section TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hr TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_rmpv, matrix_vxc, rho_ao TYPE(dft_control_type), POINTER :: dft_control @@ -2023,30 +2368,40 @@ CONTAINS CALL auxbas_pw_pool%give_back_pw(wf_r) END IF + e_density_section = cube_or_openpmd(input, str_e_density_cubes, str_e_density_openpmd, logger) + ! Write cube file with electron density - IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & - "DFT%PRINT%E_DENSITY_CUBE"), cp_p_file)) THEN + IF (e_density_section%do_output) THEN CALL section_vals_val_get(dft_section, & - keyword_name="PRINT%E_DENSITY_CUBE%DENSITY_INCLUDE", & + keyword_name=e_density_section%concat_to_relative("%DENSITY_INCLUDE"), & c_val=print_density) print_density = TRIM(print_density) - append_cube = section_get_lval(input, "DFT%PRINT%E_DENSITY_CUBE%APPEND") + ! For openPMD, refer to access modes instead of APPEND key + IF (e_density_section%grid_output == grid_output_cubes) THEN + append_cube = section_get_lval(input, e_density_section%concat_to_absolute("%APPEND")) + END IF my_pos_cube = "REWIND" IF (append_cube) THEN my_pos_cube = "APPEND" END IF ! Write the info on core densities for the interface between cp2k and the XRD code ! together with the valence density they are used to compute the form factor (Fourier transform) - xrd_interface = section_get_lval(input, "DFT%PRINT%E_DENSITY_CUBE%XRD_INTERFACE") + IF (e_density_section%grid_output == grid_output_cubes) THEN + xrd_interface = section_get_lval(input, e_density_section%concat_to_absolute("%XRD_INTERFACE")) + ELSE + ! Unimplemented for openPMD, since this does not use the regular routines + xrd_interface = .FALSE. + END IF + IF (xrd_interface) THEN !cube file only contains soft density (GAPW) IF (dft_control%qs_control%gapw) print_density = "SOFT_DENSITY" filename = "ELECTRON_DENSITY" - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & + unit_nr = cp_print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & extension=".xrd", middle_name=TRIM(filename), & file_position=my_pos_cube, log_filename=.FALSE.) - ngto = section_get_ival(input, "DFT%PRINT%E_DENSITY_CUBE%NGAUSS") + ngto = section_get_ival(input, e_density_section%concat_to_absolute("%NGAUSS")) IF (output_unit > 0) THEN INQUIRE (UNIT=unit_nr, NAME=filename) WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & @@ -2149,7 +2504,7 @@ CONTAINS DEALLOCATE (ppdens, aedens, ccdens) CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE") + e_density_section%absolute_section_key) END IF IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_DENSITY") THEN @@ -2186,10 +2541,10 @@ CONTAINS rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1) filename = "TOTAL_ELECTRON_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-total-electron-density") IF (output_unit > 0) THEN IF (.NOT. mpi_io) THEN INQUIRE (UNIT=unit_nr, NAME=filename) @@ -2197,7 +2552,7 @@ CONTAINS filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The total electron density is written in cube file format to the file:", & + "The total electron density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & "q(max) [1/Angstrom] :", q_max/angstrom, & @@ -2206,11 +2561,11 @@ CONTAINS "Total electronic charge (G-space):", rho_total, & "Total electronic charge (R-space):", rho_total_rspace END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "TOTAL ELECTRON DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL ELECTRON DENSITY", & + particles=particles, zeff=zcharge, & + stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) ! Print total spin density for spin-polarized systems IF (dft_control%nspins > 1) THEN CALL pw_zero(rho_elec_gspace) @@ -2232,18 +2587,18 @@ CONTAINS rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1) filename = "TOTAL_SPIN_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-total-spin-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The total spin density is written in cube file format to the file:", & + "The total spin density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & "q(max) [1/Angstrom] :", q_max/angstrom, & @@ -2252,11 +2607,11 @@ CONTAINS "Total spin density (G-space) :", rho_total, & "Total spin density (R-space) :", rho_total_rspace END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "TOTAL SPIN DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL SPIN DENSITY", & + particles=particles, zeff=zcharge, & + stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) END IF CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace) CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace) @@ -2273,72 +2628,71 @@ CONTAINS CALL pw_axpy(rho_r(2), rho_elec_rspace) filename = "ELECTRON_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-electron-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The sum of alpha and beta density is written in cube file format to the file:", & + "The sum of alpha and beta density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "SUM OF ALPHA AND BETA DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), & - mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SUM OF ALPHA AND BETA DENSITY", & + particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), & + mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) CALL pw_copy(rho_r(1), rho_elec_rspace) CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp) filename = "SPIN_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-spin-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The spin density is written in cube file format to the file:", & + "The spin density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "SPIN DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", & + particles=particles, zeff=zcharge, & + stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace) ELSE filename = "ELECTRON_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-electron-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The electron density is written in cube file format to the file:", & + "The electron density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) END IF - CALL cp_pw_to_cube(rho_r(1), unit_nr, "ELECTRON DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_r(1), unit_nr, "ELECTRON DENSITY", & + particles=particles, zeff=zcharge, & + stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) END IF ! nspins ELSE IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_HARD_APPROX") THEN @@ -2372,30 +2726,29 @@ CONTAINS filename = "ELECTRON_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-electron-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The electron density is written in cube file format to the file:", & + "The electron density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & "Soft electronic charge (R-space) :", rho_soft, & "Hard electronic charge (R-space) :", rho_hard, & "Total electronic charge (R-space):", rho_total_rspace END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "ELECTRON DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), & - mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "ELECTRON DENSITY", & + particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), & + mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) !------------ IF (dft_control%nspins > 1) THEN @@ -2415,29 +2768,29 @@ CONTAINS filename = "SPIN_DENSITY" mpi_io = .TRUE. - unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & - fout=mpi_filename) + unit_nr = e_density_section%print_key_unit_nr(logger, input, e_density_section%absolute_section_key, & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & + fout=mpi_filename, openpmd_basename="dft-spin-density") IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN + IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN INQUIRE (UNIT=unit_nr, NAME=filename) ELSE filename = mpi_filename END IF WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & - "The spin density is written in cube file format to the file:", & + "The spin density is written in "//e_density_section%format_name//" file format to the file:", & TRIM(filename) WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & "Soft part of the spin density :", rho_soft, & "Hard part of the spin density :", rho_hard, & "Total spin density (R-space) :", rho_total_rspace END IF - CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "SPIN DENSITY", & - particles=particles, zeff=zcharge, & - stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) - CALL cp_print_key_finished_output(unit_nr, logger, input, & - "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) + CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", & + particles=particles, zeff=zcharge, & + stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io) + CALL e_density_section%print_key_finished_output(unit_nr, logger, input, & + e_density_section%absolute_section_key, mpi_io=mpi_io) END IF ! nspins CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace) DEALLOCATE (my_Q0) diff --git a/tests/HUGE_TESTS_SUPPRESSIONS b/tests/HUGE_TESTS_SUPPRESSIONS index bc3fc2f927..8d737c4792 100644 --- a/tests/HUGE_TESTS_SUPPRESSIONS +++ b/tests/HUGE_TESTS_SUPPRESSIONS @@ -182,4 +182,7 @@ QS/regtest-fftw-wisdom/test_wisdom_import.inp QS/regtest-smeagol-2/au-1x1x4-bulk.inp QS/regtest-smeagol-2/au-h2-au-1x1x1-1V-forces.inp +# Tests produced 11.69 MiB of output. +QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp + #EOF diff --git a/tests/QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp b/tests/QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp new file mode 100644 index 0000000000..c26e2fe583 --- /dev/null +++ b/tests/QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp @@ -0,0 +1,84 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT H2O-geo-ot-pdos-lumo-comp + RUN_TYPE GEO_OPT +&END GLOBAL + +&MOTION + &GEO_OPT + MAX_ITER 1 + OPTIMIZER BFGS + &END GEO_OPT +&END MOTION + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + &END MGRID + &PRINT + &ELF_OPENPMD + ADD_LAST NUMERIC + ! Write this output with ADIOS2's BP5 engine, determined by filename + ! ending. Write the entire output into a single ADIOS2 file. + OPENPMD_EXTENSION .bp5 + STRIDE 2 2 2 + &END ELF_OPENPMD + &E_DENSITY_OPENPMD + ! Write this output with HDF5, determined by filename ending. + ! Write one separate file per Iteration, determined by expansion + ! pattern %06T. + OPENPMD_EXTENSION _%06T.h5 + STRIDE 1 1 1 + &END E_DENSITY_OPENPMD + &MO_OPENPMD + NLUMO 2 + ! Determine the extension by inline-specified TOML config. + ! Can be replaced by OPENPMD_CFG_FILE for more complex configurations. + OPENPMD_CFG "backend = 'adios2'" + OPENPMD_EXTENSION .%E + &END MO_OPENPMD + &END PRINT + &QS + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION PS + EXTRAPOLATION_ORDER 2 + &END QS + &SCF + EPS_SCF 1.0E-6 + IGNORE_CONVERGENCE_FAILURE + MAX_SCF 10 + SCF_GUESS ATOMIC + &OT + &END OT + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 + H 0.000000 -0.757136 0.520545 + H 0.000000 0.757136 0.520545 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + POTENTIAL GTH-BLYP-q6 + &END KIND + &TOPOLOGY + CONNECTIVITY GENERATE + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-openpmd/TEST_FILES.toml b/tests/QS/regtest-openpmd/TEST_FILES.toml new file mode 100644 index 0000000000..47de327722 --- /dev/null +++ b/tests/QS/regtest-openpmd/TEST_FILES.toml @@ -0,0 +1,7 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tools/do_regtest + +"H2O-geo-ot-lumo-all-openpmd.inp" = [{matcher="E_total", tol=9.0e-14, ref=-17.11058846146146}] diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 0a12ab482c..4128a21c55 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -51,6 +51,7 @@ QS/regtest-gpw-1 libint libvori !ifx QS/regtest-tddfpt-force-2 libint !ifx Fist/regtest-12 QS/regtest-ot +QS/regtest-openpmd openpmd QS/regtest-ec-force libint !ifx QS/regtest-embed libint !ifx xTB/regtest-5 diff --git a/tools/precommit/check_file_properties.py b/tools/precommit/check_file_properties.py index b88365d56b..ad23d806ba 100755 --- a/tools/precommit/check_file_properties.py +++ b/tools/precommit/check_file_properties.py @@ -76,6 +76,8 @@ FLAG_EXCEPTIONS = ( r"__LIBXSMM2", r"CPVERSION", r"_WIN32", + r"OPENPMDAPI_VERSION_GE", + r"openPMD_HAVE_MPI", # TODO: Add CMake support for the following flags or remove the corresponding code. # See also https://github.com/cp2k/cp2k/issues/4611 r"__PW_FPGA", diff --git a/tools/precommit/precommit.py b/tools/precommit/precommit.py index 7ceaebba82..845e9bef9e 100755 --- a/tools/precommit/precommit.py +++ b/tools/precommit/precommit.py @@ -225,11 +225,12 @@ def process_file(fn: str, allow_modifications: bool) -> None: run_remote_tool("clangformat", fn) if re.match(r".*\.(cc|cpp|cxx|hcc|hpp|hxx)$", fn): - if fn.endswith("/torch_c_api.cpp"): - # Begrudgingly tolerated because PyTorch has no C API. - run_remote_tool("clangformat", fn) - elif fn.endswith("/ace_c_api.cpp"): - # same as PyTorch + if ( + fn.endswith("/torch_c_api.cpp") + or fn.endswith("/ace_c_api.cpp") + or fn.endswith("/openPMD.cpp") + ): + # Begrudgingly tolerated because PyTorch/Ace/openPMD have no C API. run_remote_tool("clangformat", fn) else: raise Exception(f"C++ is not supported.") diff --git a/tools/spack/cp2k_deps_psmp.yaml b/tools/spack/cp2k_deps_psmp.yaml index f4900b2906..c6f4b3410f 100644 --- a/tools/spack/cp2k_deps_psmp.yaml +++ b/tools/spack/cp2k_deps_psmp.yaml @@ -77,6 +77,12 @@ spack: require: - +kxc - build_system=cmake + openpmd-api: + require: + - +adios2 + adios2: + require: + - ~bzip2 pexsi: require: - +fortran @@ -130,6 +136,8 @@ spack: - "libxc@7.0.0" # - "libxsmm@1.17" - "mimic-mcl@3.0.0" + - "adios2@2.10.2" + - "openpmd-api@0.16.1" - "pexsi@2.0.0" - "plumed@2.9.2" - "py-torch@2.9.0"