comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Write cube values with explicit three-digit exponents (#5245)

Browse source 
Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-19 05:38:47 +02:00 • committed by GitHub
parent 666f2d7d18
commit f658ef2d91
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 137 additions and 58 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/CMakeLists.txt
View file

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

5
src/input_cp2k_external.F
View file

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

12
src/input_cp2k_opt.F
View file

@ -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 <FILENAME>", &
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 <FILENAME>", &
type_of_var=lchar_t)
CALL section_add_keyword(section, keyword)

8
src/input_cp2k_qs.F
View file

@ -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 <FILENAME>", &
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 <FILENAME>", &
type_of_var=lchar_t)
CALL section_add_keyword(section, keyword)

118
src/pw/realspace_grid_cube.F
View file

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

46
src/pw/realspace_grid_cube_unittest.F Normal file
View file

@ -0,0 +1,46 @@
!--------------------------------------------------------------------------------------------------!
! 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 !
!--------------------------------------------------------------------------------------------------!
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

1
tests/UNIT_TESTS
View file

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

2
tools/cubecruncher/cubecruncher.f90
View file

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