From f658ef2d91a31bfaf38f742975535b3ef4a607d9 Mon Sep 17 00:00:00 2001 From: Dynamics of Condensed Matter <30792324+DCM-Uni-Paderborn@users.noreply.github.com> Date: Tue, 19 May 2026 05:38:47 +0200 Subject: [PATCH] Write cube values with explicit three-digit exponents (#5245) Co-authored-by: Thomas D. Kuehne --- src/CMakeLists.txt | 3 + src/input_cp2k_external.F | 5 +- src/input_cp2k_opt.F | 12 ++- src/input_cp2k_qs.F | 8 +- src/pw/realspace_grid_cube.F | 118 +++++++++++++++----------- src/pw/realspace_grid_cube_unittest.F | 46 ++++++++++ tests/UNIT_TESTS | 1 + tools/cubecruncher/cubecruncher.f90 | 2 +- 8 files changed, 137 insertions(+), 58 deletions(-) create mode 100644 src/pw/realspace_grid_cube_unittest.F diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8b1b8a889f..6609bb7be2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1519,6 +1519,7 @@ set(CP2K_PROGS_F common/gemm_square_unittest.F common/memory_utilities_unittest.F common/parallel_rng_types_unittest.F + pw/realspace_grid_cube_unittest.F metadyn_tools/graph.F motion/dumpdcd.F motion/xyz2dcd.F @@ -1865,6 +1866,7 @@ list( "gemm_square_unittest" "memory_utilities_unittest" "parallel_rng_types_unittest" + "realspace_grid_cube_unittest" "graph" "dumpdcd" "xyz2dcd" @@ -1880,6 +1882,7 @@ add_executable(cp2k-bin start/cp2k.F) add_executable(gemm_square_unittest common/gemm_square_unittest.F) add_executable(memory_utilities_unittest common/memory_utilities_unittest.F) add_executable(parallel_rng_types_unittest common/parallel_rng_types_unittest.F) +add_executable(realspace_grid_cube_unittest pw/realspace_grid_cube_unittest.F) add_executable(graph metadyn_tools/graph.F) add_executable(dumpdcd motion/dumpdcd.F) add_executable(xyz2dcd motion/xyz2dcd.F) diff --git a/src/input_cp2k_external.F b/src/input_cp2k_external.F index b5a4f7d8a2..b330060fef 100644 --- a/src/input_cp2k_external.F +++ b/src/input_cp2k_external.F @@ -116,7 +116,10 @@ CONTAINS !keyword for reading the external potential from cube file CALL keyword_create(keyword, __LOCATION__, name="READ_FROM_CUBE", & description="Switch for reading the external potential from file pot.cube. The values "// & - "of the potential must be on the grid points of the realspace grid.", & + "of the potential must be on the grid points of the realspace grid. "// & + "Whitespace-separated cube values are accepted. If adjacent values are written "// & + "without whitespace, each value must occupy a 13-character E13.5 field, as in "// & + "CP2K-generated cube files.", & usage="READ_FROM_CUBE T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/input_cp2k_opt.F b/src/input_cp2k_opt.F index 5624e1a7be..d7f4b8d64b 100644 --- a/src/input_cp2k_opt.F +++ b/src/input_cp2k_opt.F @@ -270,7 +270,9 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="READ_EMBED_POT_CUBE", & description="Read the embedding potential "// & - "(restart) from the cube file.", & + "(restart) from the cube file. Whitespace-separated cube values are accepted. If "// & + "adjacent values are written without whitespace, each value must occupy a "// & + "13-character E13.5 field, as in CP2K-generated cube files.", & usage="READ_EMBED_POT_CUBE .FALSE.", default_l_val=.FALSE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -285,7 +287,9 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="EMBED_CUBE_FILE_NAME", & description="Root of the file name where to read the embedding "// & - "potential (guess) as a cube.", & + "potential (guess) as a cube. Whitespace-separated cube values are accepted. If "// & + "adjacent values are written without whitespace, each value must occupy a "// & + "13-character E13.5 field, as in CP2K-generated cube files.", & usage="EMBED_CUBE_FILE_NAME ", & type_of_var=lchar_t) CALL section_add_keyword(section, keyword) @@ -293,7 +297,9 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="EMBED_SPIN_CUBE_FILE_NAME", & description="Root of the file name where to read the spin part "// & - "of the embedding potential (guess) as a cube.", & + "of the embedding potential (guess) as a cube. Whitespace-separated cube values are "// & + "accepted. If adjacent values are written without whitespace, each value must occupy "// & + "a 13-character E13.5 field, as in CP2K-generated cube files.", & usage="EMBED_SPIN_CUBE_FILE_NAME ", & type_of_var=lchar_t) CALL section_add_keyword(section, keyword) diff --git a/src/input_cp2k_qs.F b/src/input_cp2k_qs.F index e85634c924..c1940e9d65 100644 --- a/src/input_cp2k_qs.F +++ b/src/input_cp2k_qs.F @@ -512,7 +512,9 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="EMBED_CUBE_FILE_NAME", & description="Root of the file name where to read the embedding "// & - "potential (guess) as a cube.", & + "potential (guess) as a cube. Whitespace-separated cube values are accepted. If "// & + "adjacent values are written without whitespace, each value must occupy a "// & + "13-character E13.5 field, as in CP2K-generated cube files.", & usage="EMBED_CUBE_FILE_NAME ", & type_of_var=lchar_t) CALL section_add_keyword(section, keyword) @@ -520,7 +522,9 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="EMBED_SPIN_CUBE_FILE_NAME", & description="Root of the file name where to read the spin part "// & - "of the embedding potential (guess) as a cube.", & + "of the embedding potential (guess) as a cube. Whitespace-separated cube values are "// & + "accepted. If adjacent values are written without whitespace, each value must occupy "// & + "a 13-character E13.5 field, as in CP2K-generated cube files.", & usage="EMBED_SPIN_CUBE_FILE_NAME ", & type_of_var=lchar_t) CALL section_add_keyword(section, keyword) diff --git a/src/pw/realspace_grid_cube.F b/src/pw/realspace_grid_cube.F index 7d72341b1b..c470d2ac77 100644 --- a/src/pw/realspace_grid_cube.F +++ b/src/pw/realspace_grid_cube.F @@ -25,15 +25,50 @@ MODULE realspace_grid_cube PRIVATE - PUBLIC :: pw_to_cube, cube_to_pw, pw_to_simple_volumetric + PUBLIC :: cube_read_values, cube_value_format, cube_values_format, cube_to_pw, pw_to_cube, & + pw_to_simple_volumetric CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'realspace_grid_cube' + INTEGER, PARAMETER, PRIVATE :: cube_entry_len = 13, & + cube_num_entries_line = 6 + INTEGER, PARAMETER, PRIVATE :: cube_line_len = cube_entry_len*cube_num_entries_line + CHARACTER(len=*), PARAMETER :: cube_value_format = '(E13.5E3)', & + cube_values_format = '(6E13.5E3)' LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. - LOGICAL, PRIVATE :: parses_linebreaks = .FALSE., & - parse_test = .TRUE. CONTAINS +! ************************************************************************************************** +!> \brief Read cube values from a character buffer. +!> \param values value buffer from one cube line or one z-slice +!> \param buffer parsed values +! ************************************************************************************************** + SUBROUTINE cube_read_values(values, buffer) + CHARACTER(LEN=*), INTENT(IN) :: values + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: buffer + + CHARACTER(LEN=cube_entry_len) :: value + INTEGER :: i, pos, readstat + + READ (values, *, IOSTAT=readstat) buffer + IF (readstat == 0) RETURN + + pos = 1 + DO i = 1, SIZE(buffer) + IF (pos + cube_entry_len - 1 > LEN(values)) CPABORT("Unexpected end of cube data.") + value = values(pos:pos + cube_entry_len - 1) + READ (value, '(E13.5)', IOSTAT=readstat) buffer(i) + IF (readstat /= 0) CPABORT("Bad value while reading cube data.") + pos = pos + cube_entry_len + IF (MODULO(i, cube_num_entries_line) == 0) THEN + IF (pos <= LEN(values)) THEN + IF (values(pos:pos) == NEW_LINE('C')) pos = pos + 1 + END IF + END IF + END DO + + END SUBROUTINE cube_read_values + ! ************************************************************************************************** !> \brief ... !> \param pw ... @@ -215,7 +250,7 @@ CONTAINS IF (buf(I3) < 1.E-7_dp) buf(I3) = 0.0_dp END DO END IF - WRITE (unit_nr, '(6E13.5)') (buf(I3), I3=L3, U3, my_stride(3)) + WRITE (unit_nr, cube_values_format) (buf(I3), I3=L3, U3, my_stride(3)) END IF ! this double loop generates so many messages that it can overload @@ -268,9 +303,11 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'cube_to_pw' INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6 - INTEGER :: extunit, handle, i, j, k, msglen, & - my_rank, nat, ndum, num_linebreak, & - num_pe, output_unit, size_of_z, tag + CHARACTER(LEN=cube_line_len) :: value_line + INTEGER :: extunit, handle, i, j, k, last_z, & + msglen, my_rank, nat, ndum, & + num_linebreak, num_pe, output_unit, & + size_of_z, tag INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, & npoints_local, ubounds, ubounds_local LOGICAL :: be_silent @@ -344,7 +381,11 @@ CONTAINS DO i = lbounds(1), ubounds(1) DO j = lbounds(2), ubounds(2) IF (my_rank == 0) THEN - READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3)) + DO k = lbounds(3), ubounds(3), cube_num_entries_line + last_z = MIN(k + cube_num_entries_line - 1, ubounds(3)) + READ (extunit, '(A)') value_line + CALL cube_read_values(value_line, buffer(k:last_z)) + END DO END IF CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0) @@ -397,18 +438,18 @@ CONTAINS INTEGER, INTENT(in) :: msglen LOGICAL, INTENT(in), OPTIONAL :: silent - INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6 - + CHARACTER(LEN=cube_line_len) :: value_line INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, & npoints_local, ubounds, ubounds_local INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: blocklengths INTEGER(kind=file_offset), ALLOCATABLE, & DIMENSION(:), TARGET :: displacements INTEGER(kind=file_offset) :: BOF - INTEGER :: extunit_handle, i, ientry, islice, j, k, my_rank, nat, ndum, nslices, num_pe, & - offset_global, output_unit, pos, readstat, size_of_z, tag + INTEGER :: extunit_handle, i, islice, j, k, last_z, & + my_rank, nat, ndum, nslices, num_pe, & + offset_global, output_unit, size_of_z, & + tag CHARACTER(LEN=msglen), ALLOCATABLE, DIMENSION(:) :: readbuffer - CHARACTER(LEN=msglen) :: tmp LOGICAL :: be_silent, should_read(2) REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer REAL(kind=dp), DIMENSION(3) :: dr, rdum @@ -481,7 +522,9 @@ CONTAINS END IF ! Sync offset and start parallel read CALL gid%bcast(offset_global, grid%pw_grid%para%group%source) - BOF = offset_global + ! INQUIRE(POS=...) returns a 1-based stream position, whereas MPI-IO + ! file offsets are 0-based. + BOF = offset_global - 1 CALL extunit%open(groupid=gid, filepath=filename, amode_status=file_amode_rdonly) ! Determine byte offsets for each grid z-slice which are local to a process ALLOCATE (displacements(nslices)) @@ -525,39 +568,8 @@ CONTAINS j = lbounds_local(2) ALLOCATE (buffer(lbounds(3):ubounds(3))) buffer = 0.0_dp - ! Test if the compiler supports parsing lines with line breaks in them - IF (parse_test) THEN - READ (readbuffer(1), *, IOSTAT=readstat) (buffer(k), k=lbounds(3), ubounds(3)) - IF (readstat == 0) THEN - parses_linebreaks = .TRUE. - ELSE - parses_linebreaks = .FALSE. - END IF - parse_test = .FALSE. - buffer = 0.0_dp - END IF DO islice = 1, nslices - IF (parses_linebreaks) THEN - ! Use list-directed conversion if supported - ! Produces faster code, but maybe the latter method could be optimized - READ (readbuffer(islice), *) (buffer(k), k=lbounds(3), ubounds(3)) - ELSE - ! Convert values by going through the z-slice one value at a time - tmp = readbuffer(islice) - pos = 1 - ientry = 1 - DO k = lbounds_local(3), ubounds_local(3) - IF (MODULO(ientry, num_entries_line) == 0 .OR. k == ubounds_local(3)) THEN - ! Last value on line, dont read line break - READ (tmp(pos:pos + (entry_len - 2)), '(E12.5)') buffer(k) - pos = pos + (entry_len + 1) - ELSE - READ (tmp(pos:pos + (entry_len - 1)), '(E13.5)') buffer(k) - pos = pos + entry_len - END IF - ientry = ientry + 1 - END DO - END IF + CALL cube_read_values(readbuffer(islice), buffer(lbounds(3):ubounds(3))) ! Optionally scale cube file values grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3)) j = j + 1 @@ -606,7 +618,11 @@ CONTAINS DO i = lbounds(1), ubounds(1) DO j = lbounds(2), ubounds(2) IF (my_rank == 0) THEN - READ (extunit_handle, *) (buffer(k), k=lbounds(3), ubounds(3)) + DO k = lbounds(3), ubounds(3), cube_num_entries_line + last_z = MIN(k + cube_num_entries_line - 1, ubounds(3)) + READ (extunit_handle, '(A)') value_line + CALL cube_read_values(value_line, buffer(k:last_z)) + END DO END IF CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0) @@ -802,9 +818,9 @@ CONTAINS counter = 0 DO k = lbounds(3), ubounds(3), stride(3) IF (zero_tails .AND. grid%array(i, j, k) < 1.E-7_dp) THEN - WRITE (value, '(E13.5)') 0.0_dp + WRITE (value, cube_value_format) 0.0_dp ELSE - WRITE (value, '(E13.5)') grid%array(i, j, k) + WRITE (value, cube_value_format) grid%array(i, j, k) END IF tmp = TRIM(tmp)//TRIM(value) counter = counter + 1 @@ -974,7 +990,7 @@ CONTAINS pw%pw_grid%dh(3, 3)*I3 IF (unit_nr > 0) THEN - WRITE (unit_nr, '(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(I3) + WRITE (unit_nr, '(6E13.5E3, 6E13.5E3, 6E13.5E3, 6E13.5E3)') x, y, z, buf(I3) END IF END DO @@ -994,7 +1010,7 @@ CONTAINS pw%pw_grid%dh(3, 3)*I3 IF (unit_nr > 0) THEN - WRITE (unit_nr, '(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(I3), buf2(I3) + WRITE (unit_nr, '(6E13.5E3, 6E13.5E3, 6E13.5E3, 6E13.5E3)') x, y, z, buf(I3), buf2(I3) END IF END DO diff --git a/src/pw/realspace_grid_cube_unittest.F b/src/pw/realspace_grid_cube_unittest.F new file mode 100644 index 0000000000..7732cea3ea --- /dev/null +++ b/src/pw/realspace_grid_cube_unittest.F @@ -0,0 +1,46 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! +PROGRAM realspace_grid_cube_unittest + USE kinds, ONLY: dp + USE realspace_grid_cube, ONLY: cube_read_values,& + cube_value_format,& + cube_values_format + + IMPLICIT NONE + + CHARACTER(LEN=13) :: value + CHARACTER(LEN=92) :: slice_values + CHARACTER(LEN=79) :: values + REAL(KIND=dp), DIMENSION(6) :: buffer, reference + REAL(KIND=dp), DIMENSION(7) :: slice_buffer, slice_reference + + WRITE (value, cube_value_format) 0.33004E-101_dp + IF (INDEX(value, "E-101") <= 0) & + ERROR STOP "Cube value format must preserve three-digit negative exponents." + + reference = [0.0_dp, -0.91458E-35_dp, 0.33004E-101_dp, & + -0.32487E-103_dp, 1.0_dp, -1.0_dp] + WRITE (values(1:78), cube_values_format) reference + values(79:79) = NEW_LINE("C") + IF (INDEX(values, "E-101") <= 0) & + ERROR STOP "Cube line format must preserve E-101." + IF (INDEX(values, "E-103") <= 0) & + ERROR STOP "Cube line format must preserve E-103." + CALL cube_read_values(values, buffer) + IF (MAXVAL(ABS(buffer - reference)) > 1.0E-12_dp) & + ERROR STOP "Cube reader must parse adjacent explicit-exponent values." + + slice_reference(1:6) = reference + slice_reference(7) = 0.27183E-123_dp + WRITE (slice_values(1:78), cube_values_format) slice_reference(1:6) + slice_values(79:79) = NEW_LINE("C") + WRITE (slice_values(80:92), cube_value_format) slice_reference(7) + CALL cube_read_values(slice_values, slice_buffer) + IF (MAXVAL(ABS(slice_buffer - slice_reference)) > 1.0E-12_dp) & + ERROR STOP "Cube reader must parse multi-line z-slices." + +END PROGRAM realspace_grid_cube_unittest diff --git a/tests/UNIT_TESTS b/tests/UNIT_TESTS index 74063ae158..6051b88da1 100644 --- a/tests/UNIT_TESTS +++ b/tests/UNIT_TESTS @@ -9,6 +9,7 @@ gemm_square_unittest memory_utilities_unittest nequip_unittest libtorch parallel_rng_types_unittest +realspace_grid_cube_unittest gx_ac_unittest greenx #EOF diff --git a/tools/cubecruncher/cubecruncher.f90 b/tools/cubecruncher/cubecruncher.f90 index 01559b51ed..b864a7979e 100644 --- a/tools/cubecruncher/cubecruncher.f90 +++ b/tools/cubecruncher/cubecruncher.f90 @@ -113,7 +113,7 @@ CONTAINS ENDDO DO I1=1,cube%npoints(1),stride DO I2=1,cube%npoints(2),stride - WRITE(iunit,'(6E13.5)') (cube%grid(I1,I2,I3)*mult_fact,I3=1,cube%npoints(3),stride) + WRITE(iunit,'(6E13.5E3)') (cube%grid(I1,I2,I3)*mult_fact,I3=1,cube%npoints(3),stride) ENDDO ENDDO END SUBROUTINE write_cube