diff --git a/docs/advanced/elec_properties/hs_matrix.md b/docs/advanced/elec_properties/hs_matrix.md index e10c3fac95..c7285a5f96 100644 --- a/docs/advanced/elec_properties/hs_matrix.md +++ b/docs/advanced/elec_properties/hs_matrix.md @@ -18,7 +18,7 @@ Use [out_hsk](../input_files/input-main.md#out_hsk) to print the upper triangula | --- | --- | | `0` | Disabled | | `1` | Text; an optional second value controls precision, for example `out_hsk 1 12` | -| `2` | Reserved for future binary output; not implemented | +| `2` | Native binary `.dat` output | | `3` | Reserved for H(k)/S(k) NPZ output; not implemented | The legacy keyword `out_mat_hs 1 [precision]` remains supported as an alias for `out_hsk 1 [precision]`. If both names are present, `out_hsk` takes precedence. @@ -45,6 +45,18 @@ Each output block starts with a comment header containing the one-based ionic-st For multi-k calculations, the matrices are Hermitian and each matrix element is written as `(real,imag)`. For gamma-only calculations, the matrices are symmetric and the matrix elements are written as real numbers. +### Native Binary Format + +For `out_hsk 2`, the filenames in the table above use `.dat` instead of `.txt`. Each matrix record is written without padding or a self-describing header: + +1. the matrix dimension as a native C++ `int`; +2. the upper-triangular elements in row-major order, from `(0,0)` through `(0,N-1)`, then `(1,1)` through `(1,N-1)`, and so on; +3. each gamma-only element as one native `double`, or each multi-k/spinor element as two consecutive native `double` values containing the real and imaginary parts. + +The file therefore contains `sizeof(int) + N(N+1)/2 * sizeof(double)` bytes for a gamma-only record and `sizeof(int) + N(N+1) * sizeof(double)` bytes for a multi-k or spinor record. The format uses the host integer representation and byte order and is intended for readers using a compatible ABI. + +When `out_app_flag` is true, the first ionic step truncates the shared file and every later step appends another complete record. When it is false, each ionic step has a separate filename containing `g${step}` before `_nao`. + ## out_hsr The output of $H(R)$ and $S(R)$ matrices is controlled by [out_hsr](../input_files/input-main.md#out_hsr). It is available for both gamma-only and multi-k LCAO calculations: diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 0ed67ff2b7..375b451eed 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2003,7 +2003,7 @@ - **Description**: Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: - 0: disabled; - 1: text output; the optional second integer controls precision and defaults to 8; - - 2: reserved for binary output, which is not implemented yet; + - 2: binary output in the native ABACUS .dat format; - 3: NPZ output, which is not implemented for H(k)/S(k). The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. @@ -2014,6 +2014,8 @@ - Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. - Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. - Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. + For binary output, the same names use the .dat suffix. Each native binary record contains the matrix dimension as an int followed by the row-major upper triangle. Gamma-only elements are doubles; multi-k and spinor elements are pairs of doubles containing the real and imaginary parts. Native integer representation and byte order are used. + When out_app_flag is true, the first ionic step truncates the file and later steps append complete records. When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. > Note: In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 7c142a45cb..0c2494dee3 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -3112,7 +3112,7 @@ parameters: Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: * 0: disabled; * 1: text output; the optional second integer controls precision and defaults to 8; - * 2: reserved for binary output, which is not implemented yet; + * 2: binary output in the native ABACUS .dat format; * 3: NPZ output, which is not implemented for H(k)/S(k). The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. @@ -3122,6 +3122,8 @@ parameters: * Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. * Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. * Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. + For binary output, the same names use the .dat suffix. Each native binary record contains the matrix dimension as an int followed by the row-major upper triangle. Gamma-only elements are doubles; multi-k and spinor elements are pairs of doubles containing the real and imaginary parts. Native integer representation and byte order are used. + When out_app_flag is true, the first ionic step truncates the file and later steps append complete records. When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. [NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index 73ae4d6a51..cd3371f9be 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -199,7 +199,8 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, //------------------------------------------------------------------ // 4) Output H(k) and S(k) matrices for each k-point //------------------------------------------------------------------ - if (inp.out_hsk[0] == 1) + const int hsk_out_type = inp.out_hsk[0]; + if (hsk_out_type == 1 || hsk_out_type == 2) { const int precision = inp.out_hsk[1]; ModuleIO::write_hsk(global_out_dir, @@ -213,6 +214,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, gamma_only, out_app_flag, istep, + hsk_out_type, precision, GlobalV::ofs_running); } diff --git a/source/source_io/module_hs/write_HS.h b/source/source_io/module_hs/write_HS.h index f2b0d05b8a..c511a2259a 100644 --- a/source/source_io/module_hs/write_HS.h +++ b/source/source_io/module_hs/write_HS.h @@ -26,6 +26,7 @@ namespace ModuleIO const bool gamma_only, const bool out_app_flag, const int istep, + const int out_type, const int precision, std::ofstream &ofs_running); diff --git a/source/source_io/module_hs/write_HS.hpp b/source/source_io/module_hs/write_HS.hpp index 86edca54e1..c3bd846a55 100644 --- a/source/source_io/module_hs/write_HS.hpp +++ b/source/source_io/module_hs/write_HS.hpp @@ -21,6 +21,7 @@ void ModuleIO::write_hsk( const bool gamma_only, const bool out_app_flag, const int istep, + const int out_type, const int precision, std::ofstream &ofs_running) { @@ -40,25 +41,21 @@ void ModuleIO::write_hsk( for (int ik = 0; ik < nks; ++ik) { p_hamilt->updateHk(ik); - bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core - // calculation, noted by zhengdy-soc + const bool binary = (out_type == 2); hamilt::MatrixBlock h_mat; hamilt::MatrixBlock s_mat; p_hamilt->matrix(h_mat, s_mat); - const int out_label=1; // 1: .txt, 2: .dat - std::string h_fn = ModuleIO::filename_output(global_out_dir, "hk","nao",ik,ik2iktot,nspin,nkstot, - out_label,out_app_flag,gamma_only,istep); + out_type,out_app_flag,gamma_only,istep); ModuleIO::save_mat(istep, h_mat.p, PARAM.globalv.nlocal, - bit, + binary, precision, 1, out_app_flag, @@ -77,14 +74,14 @@ void ModuleIO::write_hsk( std::string s_fn = ModuleIO::filename_output(global_out_dir, "sk","nao",ik,ik2iktot,nspin,nkstot, - out_label,out_app_flag,gamma_only,istep); + out_type,out_app_flag,gamma_only,istep); ofs_running << " The output filename is " << s_fn << std::endl; ModuleIO::save_mat(istep, s_mat.p, PARAM.globalv.nlocal, - bit, + binary, precision, 1, out_app_flag, @@ -123,7 +120,8 @@ void ModuleIO::save_mat(const int istep, if (drank == 0) { - out_matrix = fopen(filename.c_str(), "wb"); + const char* mode = (app && istep > 0) ? "ab" : "wb"; + out_matrix = fopen(filename.c_str(), mode); if (out_matrix == nullptr) { ModuleBase::WARNING_QUIT("ModuleIO::save_mat", "Cannot open matrix file: " + filename); @@ -184,7 +182,8 @@ void ModuleIO::save_mat(const int istep, } // write .dat file without MPI #else - FILE* out_matrix = fopen(filename.c_str(), "wb"); + const char* mode = (app && istep > 0) ? "ab" : "wb"; + FILE* out_matrix = fopen(filename.c_str(), mode); if (out_matrix == nullptr) { ModuleBase::WARNING_QUIT("ModuleIO::save_mat", "Cannot open matrix file: " + filename); diff --git a/source/source_io/module_parameter/read_input_item_output.cpp b/source/source_io/module_parameter/read_input_item_output.cpp index 2818b9b523..15868d0eb4 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -534,7 +534,7 @@ Also controled by out_freq_ion and out_app_flag. item.description = R"(Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: * 0: disabled; * 1: text output; the optional second integer controls precision and defaults to 8; -* 2: reserved for binary output, which is not implemented yet; +* 2: binary output in the native ABACUS .dat format; * 3: NPZ output, which is not implemented for H(k)/S(k). The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. @@ -544,6 +544,8 @@ The output is also controlled by out_freq_ion and out_app_flag. For more informa * Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. * Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. * Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. +For binary output, the same names use the .dat suffix. Each native binary record contains the matrix dimension as an int followed by the row-major upper triangle. Gamma-only elements are doubles; multi-k and spinor elements are pairs of doubles containing the real and imaginary parts. Native integer representation and byte order are used. +When out_app_flag is true, the first ionic step truncates the file and later steps append complete records. When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. [NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc.)"; @@ -576,10 +578,6 @@ When out_app_flag is false, g followed by the one-based ionic-step index is inse { ModuleBase::WARNING_QUIT("ReadInput", "out_hsk format must be 0, 1, 2, or 3"); } - if (format == 2) - { - ModuleBase::WARNING_QUIT("ReadInput", "out_hsk binary output is reserved but not implemented"); - } if (format == 3) { ModuleBase::WARNING_QUIT("ReadInput", "out_hsk NPZ output is not implemented"); diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index e0c84b8468..04903ec1a5 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -53,6 +53,18 @@ AddTest( SOURCES binstream_test.cpp ../module_output/binstream.cpp ) +AddTest( + TARGET MODULE_IO_write_hsk_binary_test + LIBS parameter base device + SOURCES write_hsk_binary_test.cpp ../../source_basis/module_ao/parallel_orbitals.cpp +) + +add_test( + NAME MODULE_IO_write_hsk_binary_test_parallel + COMMAND mpirun -np 2 ./MODULE_IO_write_hsk_binary_test + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) + AddTest( TARGET MODULE_IO_write_eig_occ_test LIBS parameter base device symmetry diff --git a/source/source_io/test/write_hsk_binary_test.cpp b/source/source_io/test/write_hsk_binary_test.cpp new file mode 100644 index 0000000000..ff16a7aa7f --- /dev/null +++ b/source/source_io/test/write_hsk_binary_test.cpp @@ -0,0 +1,215 @@ +#include "source_io/module_hs/write_HS.h" + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include + +#ifdef __MPI +#include +#endif + +namespace +{ + +int test_rank = 0; + +void initialize_distribution(Parallel_Orbitals& pv, const int dim) +{ +#ifdef __MPI + ASSERT_EQ(pv.init(dim, dim, 1, DIAG_WORLD), 0); +#else + pv.set_serial(dim, dim); +#endif +} + +template +std::vector distribute_matrix(const Parallel_Orbitals& pv, const std::vector& global, const int dim) +{ + std::vector local(pv.get_local_size(), T()); + for (int i = 0; i < dim; ++i) + { + const int ir = pv.global2local_row(i); + if (ir < 0) + { + continue; + } + for (int j = 0; j < dim; ++j) + { + const int ic = pv.global2local_col(j); + if (ic >= 0) + { + local[ir * pv.ncol + ic] = global[i * dim + j]; + } + } + } + return local; +} + +template +std::vector upper_triangle(const std::vector& matrix, const int dim) +{ + std::vector values; + values.reserve(dim * (dim + 1) / 2); + for (int i = 0; i < dim; ++i) + { + for (int j = i; j < dim; ++j) + { + values.push_back(matrix[i * dim + j]); + } + } + return values; +} + +template +std::vector read_record(std::ifstream& ifs, const int expected_dim) +{ + int dim = 0; + ifs.read(reinterpret_cast(&dim), sizeof(int)); + EXPECT_TRUE(ifs.good()); + EXPECT_EQ(dim, expected_dim); + + std::vector values(expected_dim * (expected_dim + 1) / 2); + ifs.read(reinterpret_cast(values.data()), values.size() * sizeof(T)); + EXPECT_TRUE(ifs.good()); + return values; +} + +void remove_test_file(const std::string& filename) +{ + if (test_rank == 0) + { + std::remove(filename.c_str()); + } +#ifdef __MPI + MPI_Barrier(DIAG_WORLD); +#endif +} + +} // namespace + +TEST(WriteHskBinary, WritesGammaUpperTriangle) +{ + const int dim = 3; + const std::string filename = "write_hsk_gamma.dat"; + remove_test_file(filename); + + Parallel_Orbitals pv; + initialize_distribution(pv, dim); + std::vector global(dim * dim); + for (int i = 0; i < dim * dim; ++i) + { + global[i] = i + 0.5; + } + const std::vector local = distribute_matrix(pv, global, dim); + + ModuleIO::save_mat(0, local.data(), dim, true, 8, true, false, filename, pv, test_rank); + + if (test_rank == 0) + { + std::ifstream ifs(filename.c_str(), std::ios::binary | std::ios::ate); + ASSERT_TRUE(ifs.is_open()); + EXPECT_EQ(ifs.tellg(), static_cast(sizeof(int) + 6 * sizeof(double))); + ifs.seekg(0); + EXPECT_EQ(read_record(ifs, dim), upper_triangle(global, dim)); + EXPECT_EQ(ifs.peek(), std::ifstream::traits_type::eof()); + } + + remove_test_file(filename); +} + +TEST(WriteHskBinary, WritesComplexUpperTriangleWithMpiReduction) +{ + const int dim = 4; + const std::string filename = "write_hsk_complex.dat"; + remove_test_file(filename); + + Parallel_Orbitals pv; + initialize_distribution(pv, dim); + std::vector> global(dim * dim); + for (int i = 0; i < dim * dim; ++i) + { + global[i] = std::complex(i + 0.25, -i - 0.75); + } + const std::vector> local = distribute_matrix(pv, global, dim); + + ModuleIO::save_mat(0, local.data(), dim, true, 8, true, false, filename, pv, test_rank); + + if (test_rank == 0) + { + std::ifstream ifs(filename.c_str(), std::ios::binary | std::ios::ate); + ASSERT_TRUE(ifs.is_open()); + const int element_count = dim * (dim + 1) / 2; + EXPECT_EQ(ifs.tellg(), static_cast(sizeof(int) + element_count * sizeof(std::complex))); + ifs.seekg(0); + EXPECT_EQ(read_record>(ifs, dim), upper_triangle(global, dim)); + EXPECT_EQ(ifs.peek(), std::ifstream::traits_type::eof()); + } + + remove_test_file(filename); +} + +TEST(WriteHskBinary, AppendsCompleteRecordsAndCanOverwrite) +{ + const int dim = 2; + const std::string filename = "write_hsk_append.dat"; + remove_test_file(filename); + + Parallel_Orbitals pv; + initialize_distribution(pv, dim); + const std::vector first = {1.0, 2.0, 3.0, 4.0}; + const std::vector second = {5.0, 6.0, 7.0, 8.0}; + const std::vector replacement = {9.0, 10.0, 11.0, 12.0}; + const std::vector first_local = distribute_matrix(pv, first, dim); + const std::vector second_local = distribute_matrix(pv, second, dim); + const std::vector replacement_local = distribute_matrix(pv, replacement, dim); + + ModuleIO::save_mat(0, first_local.data(), dim, true, 8, true, true, filename, pv, test_rank); + ModuleIO::save_mat(1, second_local.data(), dim, true, 8, true, true, filename, pv, test_rank); + + if (test_rank == 0) + { + std::ifstream ifs(filename.c_str(), std::ios::binary); + ASSERT_TRUE(ifs.is_open()); + EXPECT_EQ(read_record(ifs, dim), upper_triangle(first, dim)); + EXPECT_EQ(read_record(ifs, dim), upper_triangle(second, dim)); + EXPECT_EQ(ifs.peek(), std::ifstream::traits_type::eof()); + } + +#ifdef __MPI + MPI_Barrier(DIAG_WORLD); +#endif + ModuleIO::save_mat(2, replacement_local.data(), dim, true, 8, true, false, filename, pv, test_rank); + + if (test_rank == 0) + { + std::ifstream ifs(filename.c_str(), std::ios::binary); + ASSERT_TRUE(ifs.is_open()); + EXPECT_EQ(read_record(ifs, dim), upper_triangle(replacement, dim)); + EXPECT_EQ(ifs.peek(), std::ifstream::traits_type::eof()); + } + + remove_test_file(filename); +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); + DIAG_WORLD = MPI_COMM_WORLD; + MPI_Comm_rank(DIAG_WORLD, &test_rank); +#else + test_rank = 0; +#endif + + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index c59e96e484..7bde1de9b7 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -1026,11 +1026,11 @@ TEST_F(InputTest, Item_test) EXPECT_EQ(param.input.out_hsk[0], 1); EXPECT_EQ(param.input.out_hsk[1], 12); - param.input.out_hsk[0] = 2; - testing::internal::CaptureStdout(); - EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("reserved but not implemented")); + it->second.str_values = {"2", "12"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_hsk[0], 2); + EXPECT_EQ(param.input.out_hsk[1], 12); + it->second.check_value(it->second, param); param.input.out_hsk[0] = 3; testing::internal::CaptureStdout(); diff --git a/tests/03_NAO_multik/CASES_CPU.txt b/tests/03_NAO_multik/CASES_CPU.txt index beb5b15163..11193e32b9 100644 --- a/tests/03_NAO_multik/CASES_CPU.txt +++ b/tests/03_NAO_multik/CASES_CPU.txt @@ -30,6 +30,7 @@ scf_out_chg_tau scf_out_pot3 scf_out_dmr_dmk scf_out_hsk +scf_out_hsk_binary scf_out_hsr scf_out_hsr_spin4 scf_out_dh_t diff --git a/tests/03_NAO_multik/scf_out_hsk_binary/INPUT b/tests/03_NAO_multik/scf_out_hsk_binary/INPUT new file mode 100644 index 0000000000..c18698d7bf --- /dev/null +++ b/tests/03_NAO_multik/scf_out_hsk_binary/INPUT @@ -0,0 +1,29 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation scf + +nbands 6 +symmetry 0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +gamma_only 0 + +#Parameters (2.Iteration) +ecutwfc 20 +scf_thr 1e-8 +scf_nmax 50 + +#Parameters (3.Basis) +basis_type lcao + +#Parameters (4.Smearing) +smearing_method gauss +smearing_sigma 0.002 + +#Parameters (5.Mixing) +mixing_type broyden +mixing_beta 0.7 + +out_hsk 2 +ks_solver scalapack_gvx diff --git a/tests/03_NAO_multik/scf_out_hsk_binary/KPT b/tests/03_NAO_multik/scf_out_hsk_binary/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/03_NAO_multik/scf_out_hsk_binary/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/03_NAO_multik/scf_out_hsk_binary/README b/tests/03_NAO_multik/scf_out_hsk_binary/README new file mode 100644 index 0000000000..377ca172cc --- /dev/null +++ b/tests/03_NAO_multik/scf_out_hsk_binary/README @@ -0,0 +1 @@ +Test native binary H(k)/S(k) output against the existing text matrix references. diff --git a/tests/03_NAO_multik/scf_out_hsk_binary/STRU b/tests/03_NAO_multik/scf_out_hsk_binary/STRU new file mode 100644 index 0000000000..6a4ecd0ee8 --- /dev/null +++ b/tests/03_NAO_multik/scf_out_hsk_binary/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Si 14 Si_dojo_nsoc.upf + +NUMERICAL_ORBITAL +Si_dojo_6au_sz.orb + +LATTICE_CONSTANT +10.2 // add lattice constant + +LATTICE_VECTORS +0.0 0.5 0.5 +0.5 0.0 0.5 +0.5 0.5 0.0 + +ATOMIC_POSITIONS +Direct + +Si // Element type +0.0 // magnetism +2 +0.00 0.00 0.00 1 1 1 +0.25 0.25 0.25 1 1 1 diff --git a/tests/03_NAO_multik/scf_out_hsk_binary/result.ref b/tests/03_NAO_multik/scf_out_hsk_binary/result.ref new file mode 100644 index 0000000000..e92647c437 --- /dev/null +++ b/tests/03_NAO_multik/scf_out_hsk_binary/result.ref @@ -0,0 +1,5 @@ +etotref -194.1602102900063 +etotperatomref -97.0801051450 +CompareH_pass 0 +CompareS_pass 0 +totaltimeref 0.33 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index e860f1fae4..d9591d882b 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -372,6 +372,21 @@ if ! test -z "$has_hs" && [ $has_hs == 1 ]; then echo "CompareS_pass $?" >>$1 fi fi +elif ! test -z "$has_hs" && [ $has_hs == 2 ]; then + HSK_BINARY_COMPARE="../../integrate/tools/compare_hsk_binary.py" + if ! test -z "$gamma_only" && [ $gamma_only == 1 ]; then + HSK_TEXT_REFERENCE_DIR="../scf_out_hk" + python3 $HSK_BINARY_COMPARE OUT.autotest/hk_nao.dat "$HSK_TEXT_REFERENCE_DIR/hk_nao.txt.ref" real 3 + echo "CompareH_pass $?" >>$1 + python3 $HSK_BINARY_COMPARE OUT.autotest/sk_nao.dat "$HSK_TEXT_REFERENCE_DIR/sk_nao.txt.ref" real 3 + echo "CompareS_pass $?" >>$1 + else + HSK_TEXT_REFERENCE_DIR="../scf_out_hsk" + python3 $HSK_BINARY_COMPARE OUT.autotest/hk2_nao.dat "$HSK_TEXT_REFERENCE_DIR/hk2_nao.txt.ref" complex 3 + echo "CompareH_pass $?" >>$1 + python3 $HSK_BINARY_COMPARE OUT.autotest/sk2_nao.dat "$HSK_TEXT_REFERENCE_DIR/sk2_nao.txt.ref" complex 3 + echo "CompareS_pass $?" >>$1 + fi fi #-------------------------------- diff --git a/tests/integrate/tools/compare_hsk_binary.py b/tests/integrate/tools/compare_hsk_binary.py new file mode 100644 index 0000000000..69d78e32a0 --- /dev/null +++ b/tests/integrate/tools/compare_hsk_binary.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 + +import re +import struct +import sys + + +def read_reference(filename, value_type): + values = [] + complex_pattern = re.compile(r"\(([^,]+),([^\)]+)\)") + with open(filename, "r", encoding="utf-8") as stream: + for line in stream: + stripped = line.strip() + if not stripped or stripped.startswith("#") or stripped.startswith("Row "): + continue + if value_type == "complex": + values.extend(complex(float(real), float(imag)) for real, imag in complex_pattern.findall(stripped)) + else: + values.extend(float(token) for token in stripped.split()) + return values + + +def read_binary(filename, value_type): + with open(filename, "rb") as stream: + header = stream.read(struct.calcsize("@i")) + if len(header) != struct.calcsize("@i"): + raise ValueError("missing native int matrix dimension") + dim = struct.unpack("@i", header)[0] + count = dim * (dim + 1) // 2 + element_format = "@dd" if value_type == "complex" else "@d" + element_size = struct.calcsize(element_format) + payload = stream.read(count * element_size) + if len(payload) != count * element_size: + raise ValueError("binary payload is shorter than the declared upper triangle") + if stream.read(1): + raise ValueError("unexpected trailing data after the first matrix record") + + if value_type == "complex": + return [complex(real, imag) for real, imag in struct.iter_unpack(element_format, payload)] + return [value[0] for value in struct.iter_unpack(element_format, payload)] + + +def main(): + if len(sys.argv) != 5 or sys.argv[3] not in ("real", "complex"): + print("usage: compare_hsk_binary.py BINARY REFERENCE real|complex ACCURACY") + return 2 + + binary_filename, reference_filename, value_type = sys.argv[1:4] + tolerance = 10.0 ** (-int(sys.argv[4])) + try: + actual = read_binary(binary_filename, value_type) + reference = read_reference(reference_filename, value_type) + except (OSError, ValueError, struct.error) as error: + print("failed to read H(k)/S(k) output: {}".format(error)) + return 1 + + if len(actual) != len(reference): + print("matrix element count differs: binary={} reference={}".format(len(actual), len(reference))) + return 1 + + for index, (actual_value, reference_value) in enumerate(zip(actual, reference)): + if abs(actual_value - reference_value) > tolerance: + print( + "matrix element {} differs: binary={} reference={} tolerance={}".format( + index, actual_value, reference_value, tolerance + ) + ) + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main())