comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 13:45:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add openPMD output

Browse source
This commit is contained in:
Franz Pöschel 2026-01-06 12:57:32 +01:00 • committed by GitHub
parent 731b0ae562
commit 7c648bcc1f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
21 changed files with 4258 additions and 245 deletions
Download patch file Download diff file Expand all files Collapse all files

19
CMakeLists.txt
View file

@ -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()

1
cmake/cmake_cp2k.sh
View file

@ -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=$?

12
docs/technologies/libraries.md
View file

@ -229,6 +229,18 @@ greenX - Open-source file format and library. Support for greenX can be enabled
- For more information see <https://github.com/tblite/tblite>
- 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 <https://github.com/openPMD/openPMD-api/>, a equal to or
greater than 0.16.1 is required.
- For more information see <https://openpmd-api.readthedocs.io>.
- 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.

17
src/CMakeLists.txt
View file

@ -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.
$<$<BOOL:${CP2K_USE_LIBTORCH}>:${TORCH_LIBRARIES}>
$<$<BOOL:${CP2K_USE_OPENPMD}>:openPMD::openPMD>
$<$<AND:$<CONFIG:COVERAGE>,$<Fortran_COMPILER_ID:GNU>>:gcov>)
string(TIMESTAMP CP2K_TIMESTAMP "%Y-%m-%d %H:%M:%S")
@ -1790,6 +1792,7 @@ target_compile_definitions(
$<$<BOOL:${CP2K_USE_DLAF}>:__DLAF>
$<$<BOOL:${CP2K_USE_LIBXC}>:__LIBXC>
$<$<BOOL:${CP2K_USE_HDF5}>:__HDF5>
$<$<BOOL:${CP2K_USE_OPENPMD}>:__OPENPMD>
$<$<BOOL:${CP2K_USE_TREXIO}>:__TREXIO>
$<$<BOOL:${CP2K_USE_GREENX}>:__GREENX>
$<$<OR:$<BOOL:${CP2K_USE_FFTW3_}>,$<BOOL:${CP2K_USE_FFTW3_MKL_}>>:__FFTW3>

4
src/cp2k_info.F
View file

@ -274,6 +274,10 @@ CONTAINS
flags = TRIM(flags)//" greenx"
#endif
#if defined(__OPENPMD)
flags = TRIM(flags)//" openpmd"
#endif
END FUNCTION cp2k_flags
! **************************************************************************************************

90
src/cp_realspace_grid_openpmd.F Normal file
View file

@ -0,0 +1,90 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! 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

2
src/f77_interface.F
View file

@ -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

740
src/input/cp_output_handling_openpmd.F Normal file
View file

@ -0,0 +1,740 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! 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

273
src/input_cp2k_print_dft.F
View file

@ -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 ...

9
src/mpiwrap/message_passing.F
View file

@ -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)

1044
src/openPMD.F Normal file
View file

File diff suppressed because it is too large Load diff

1003
src/openPMD.cpp Normal file
View file

File diff suppressed because it is too large Load diff

526
src/pw/realspace_grid_openpmd.F Normal file
View file

@ -0,0 +1,526 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! 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

647
src/qs_scf_post_gpw.F
View file

File diff suppressed because it is too large Load diff

3
tests/HUGE_TESTS_SUPPRESSIONS
View file

@ -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

84
tests/QS/regtest-openpmd/H2O-geo-ot-lumo-all-openpmd.inp Normal file
View file

@ -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

7
tests/QS/regtest-openpmd/TEST_FILES.toml Normal file
View file

@ -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}]

1
tests/TEST_DIRS
View file

@ -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

2
tools/precommit/check_file_properties.py
View file

@ -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",

11
tools/precommit/precommit.py
View file

@ -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.")

8
tools/spack/cp2k_deps_psmp.yaml
View file

@ -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"