diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index a9eb2223762..f98dc1dc474 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -180,8 +180,13 @@ ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator( const int rank) { - // ofs_running this->ofs_ = ptr_log; + if (this->ofs_ == nullptr) + { + this->fallback_ofs_.open("/dev/null"); + this->ofs_ = &this->fallback_ofs_; + } + *ofs_ << "\n\n\n\n"; *ofs_ << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; *ofs_ << " | |" << std::endl; @@ -259,7 +264,7 @@ void ModuleIO::AngularMomentumCalculator::kernel( const char dir, const int precision) { - if (!ofs->is_open()) + if (ofs == nullptr || !ofs->is_open()) { return; } @@ -381,4 +386,4 @@ void ModuleIO::AngularMomentumCalculator::calculate( this->kernel(&ofout, ucell, d, precision); ofout.close(); } -} \ No newline at end of file +} diff --git a/source/source_io/module_hs/cal_pLpR.h b/source/source_io/module_hs/cal_pLpR.h index 9ecd15a6120..c4694596cec 100644 --- a/source/source_io/module_hs/cal_pLpR.h +++ b/source/source_io/module_hs/cal_pLpR.h @@ -88,6 +88,7 @@ #include #include #include +#include #include #include "source_cell/unitcell.h" #include "source_basis/module_nao/two_center_integrator.h" @@ -218,8 +219,9 @@ namespace ModuleIO const int rank = 0); private: + std::ofstream fallback_ofs_; // ofsrunning - std::ofstream* ofs_; + std::ofstream* ofs_ = nullptr; // the two-center-integrator std::unique_ptr calculator_; // the spherical bessel transformer @@ -246,4 +248,4 @@ namespace ModuleIO const char dir = 'x', const int precision = 10); }; -} // namespace ModuleIO \ No newline at end of file +} // namespace ModuleIO diff --git a/source/source_io/module_hs/write_HS.hpp b/source/source_io/module_hs/write_HS.hpp index 0d4950ef9c9..58e18a07554 100644 --- a/source/source_io/module_hs/write_HS.hpp +++ b/source/source_io/module_hs/write_HS.hpp @@ -280,11 +280,11 @@ void ModuleIO::save_mat(const int istep, #else if (app) { - std::ofstream out_matrix(filename.c_str(), std::ofstream::app); + out_matrix.open(filename.c_str(), std::ofstream::app); } else { - std::ofstream out_matrix(filename.c_str()); + out_matrix.open(filename.c_str()); } out_matrix << dim; diff --git a/source/source_io/module_hs/write_HS_R.cpp b/source/source_io/module_hs/write_HS_R.cpp index 2a16db86d1f..c8b0cc65ed1 100644 --- a/source/source_io/module_hs/write_HS_R.cpp +++ b/source/source_io/module_hs/write_HS_R.cpp @@ -8,34 +8,6 @@ #include "source_lcao/spar_st.h" #include "write_HS_sparse.h" -namespace { -// Helper: Convert sparse map to HContainer -template -hamilt::HContainer* sparse_map_to_hcontainer( - const std::map, std::map>>& sparse_map, - const Parallel_Orbitals& pv, - const int nbasis) -{ - hamilt::HContainer* hc = new hamilt::HContainer(&pv); - hc->set_zero(); - - for (const auto& r_entry : sparse_map) - { - const auto& R = r_entry.first; - for (const auto& row_entry : r_entry.second) - { - const size_t row = row_entry.first; - for (const auto& col_entry : row_entry.second) - { - hc->set_value(R.x, R.y, R.z, row, col_entry.first, col_entry.second); - } - } - } - - return hc; -} -} // anonymous namespace - // if 'binary=true', output binary file. // The 'sparse_thr' is the accuracy of the sparse matrix. // If the absolute value of the matrix element is less than or equal to the diff --git a/source/source_io/module_hs/write_HS_sparse.cpp b/source/source_io/module_hs/write_HS_sparse.cpp index adcdd74951d..73572aead3b 100644 --- a/source/source_io/module_hs/write_HS_sparse.cpp +++ b/source/source_io/module_hs/write_HS_sparse.cpp @@ -6,6 +6,7 @@ #include "source_lcao/module_rt/td_info.h" #include "single_R_io.h" +#include void ModuleIO::save_dH_sparse(const int& istep, const Parallel_Orbitals& pv, @@ -421,16 +422,17 @@ void ModuleIO::save_sparse( std::ofstream ofs; if (!reduce || GlobalV::DRANK == 0) { if (binary) { - int nlocal = PARAM.globalv.nlocal; + const int step = std::max(istep, 0); + const int nlocal = PARAM.globalv.nlocal; if (PARAM.inp.calculation == "md" && PARAM.inp.out_app_flag && istep) { ofs.open(sss.str().c_str(), std::ios::binary | std::ios::app); } else { ofs.open(sss.str().c_str(), std::ios::binary); } - ofs.write(reinterpret_cast(0), sizeof(int)); - ofs.write(reinterpret_cast(&nlocal), sizeof(int)); - ofs.write(reinterpret_cast(&output_R_number), sizeof(int)); + ofs.write(reinterpret_cast(&step), sizeof(int)); + ofs.write(reinterpret_cast(&nlocal), sizeof(int)); + ofs.write(reinterpret_cast(&output_R_number), sizeof(int)); } else { if (PARAM.inp.calculation == "md" && PARAM.inp.out_app_flag && istep) { diff --git a/source/source_io/module_hs/write_vxc_r.hpp b/source/source_io/module_hs/write_vxc_r.hpp index 4fe994d44a2..f5c5a320b2c 100644 --- a/source/source_io/module_hs/write_vxc_r.hpp +++ b/source/source_io/module_hs/write_vxc_r.hpp @@ -57,11 +57,10 @@ void write_Vxc_R(const int nspin, double vtxc = 0.0; // elecstate::PotXC* potxc(&rho_basis, &etxc, vtxc, nullptr); // potxc.cal_v_eff(&chg, &ucell, vr_xc); - elecstate::Potential* potxc - = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc); + elecstate::Potential potxc(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &solvent, &etxc, &vtxc); std::vector compnents_list = {"xc"}; - potxc->pot_register(compnents_list); - potxc->update_from_charge(&chg, &ucell); + potxc.pot_register(compnents_list); + potxc.update_from_charge(&chg, &ucell); // 2. allocate H(R) // (the number of hR: 1 for nspin=1, 4; 2 for nspin=2) @@ -89,19 +88,18 @@ void write_Vxc_R(const int nspin, // 3. calculate the Vxc(R) hamilt::HS_Matrix_K vxc_k_ao(pv, 1); // only hk is needed, sk is skipped - std::vector>*> vxcs_op_ao(nspin0); for (int is = 0; is < nspin0; ++is) { - vxcs_op_ao[is] = new hamilt::Veff>(&vxc_k_ao, - kv.kvec_d, - potxc, - &vxcs_R_ao[is], - &ucell, - orb_cutoff, - &gd, - nspin); - vxcs_op_ao[is]->set_current_spin(is); - vxcs_op_ao[is]->contributeHR(); + hamilt::Veff> vxcs_op_ao(&vxc_k_ao, + kv.kvec_d, + &potxc, + &vxcs_R_ao[is], + &ucell, + orb_cutoff, + &gd, + nspin); + vxcs_op_ao.set_current_spin(is); + vxcs_op_ao.contributeHR(); #ifdef __EXX if (GlobalC::exx_info.info_global.cal_exx) { diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index e9afa10db23..6dd23247472 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -283,6 +283,23 @@ AddTest( ../../source_basis/module_nao/real_gaunt_table.cpp ../../source_basis/module_nao/radial_collection.cpp ) + +AddTest( + TARGET MODULE_IO_write_hs_r_compat_test + LIBS parameter base ${math_libs} device hcontainer + SOURCES + write_hs_r_compat_test.cpp + ../module_hs/write_HS_R.cpp + ../module_hs/write_HS_sparse.cpp + ../module_hs/single_R_io.cpp + ../module_dm/write_dmr.cpp + ../module_output/ucell_io.cpp + ../module_output/sparse_matrix.cpp + ../module_output/csr_reader.cpp + ../module_output/file_reader.cpp + ../../source_basis/module_ao/parallel_orbitals.cpp + ../../source_lcao/module_hcontainer/test/tmp_mocks.cpp +) endif() AddTest( diff --git a/source/source_io/test/write_hs_r_compat_test.cpp b/source/source_io/test/write_hs_r_compat_test.cpp new file mode 100644 index 00000000000..da0bf8cc976 --- /dev/null +++ b/source/source_io/test/write_hs_r_compat_test.cpp @@ -0,0 +1,335 @@ +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#define private public +#include "source_io/module_parameter/parameter.h" +#undef private + +#include "source_base/global_variable.h" +#include "source_cell/module_neighbor/sltk_grid_driver.h" +#include "source_io/module_dm/write_dmr.h" +#include "source_io/module_hs/write_HS_R.h" +#include "source_io/module_hs/write_HS_sparse.h" +#include "source_lcao/module_hcontainer/atom_pair.h" +#include "source_lcao/module_hcontainer/hcontainer.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sparse_format +{ +void cal_dH(const UnitCell&, + const Parallel_Orbitals&, + LCAO_HS_Arrays&, + const Grid_Driver&, + const TwoCenterBundle&, + const LCAO_Orbitals&, + const int&, + const double&, + const ModuleBase::matrix&) +{ + FAIL() << "cal_dH should not be called by writer compatibility tests."; +} + +void cal_dS(const UnitCell&, + const Parallel_Orbitals&, + LCAO_HS_Arrays&, + const Grid_Driver&, + const TwoCenterBundle&, + const LCAO_Orbitals&, + const double&) +{ + FAIL() << "cal_dS should not be called by writer compatibility tests."; +} + +void cal_TR(const UnitCell&, + const Parallel_Orbitals&, + LCAO_HS_Arrays&, + const Grid_Driver&, + const TwoCenterBundle&, + const LCAO_Orbitals&, + const double&) +{ + FAIL() << "cal_TR should not be called by writer compatibility tests."; +} + +template +void cal_SR(const Parallel_Orbitals&, + std::set>&, + std::map, std::map>>&, + std::map, std::map>>>&, + const Grid_Driver&, + const double&, + hamilt::Hamilt*) +{ + FAIL() << "cal_SR should not be called by writer compatibility tests."; +} + +void destroy_dH_R_sparse(LCAO_HS_Arrays&) {} +void destroy_HS_R_sparse(LCAO_HS_Arrays&) {} +void destroy_T_R_sparse(LCAO_HS_Arrays&) {} + +template void cal_SR( + const Parallel_Orbitals&, + std::set>&, + std::map, std::map>>&, + std::map, std::map>>>&, + const Grid_Driver&, + const double&, + hamilt::Hamilt*); + +template void cal_SR>( + const Parallel_Orbitals&, + std::set>&, + std::map, std::map>>&, + std::map, std::map>>>&, + const Grid_Driver&, + const double&, + hamilt::Hamilt>*); +} // namespace sparse_format + +namespace +{ +std::string read_file(const std::string& filename) +{ + std::ifstream ifs(filename.c_str()); + std::ostringstream oss; + oss << ifs.rdbuf(); + return oss.str(); +} + +std::vector read_binary_ints(const std::string& filename, const size_t count) +{ + std::ifstream ifs(filename.c_str(), std::ios::binary); + std::vector values(count, 0); + for (size_t i = 0; i < count; ++i) + { + ifs.read(reinterpret_cast(&values[i]), sizeof(int)); + } + return values; +} + +int count_substr(const std::string& text, const std::string& pattern) +{ + int count = 0; + std::string::size_type pos = 0; + while ((pos = text.find(pattern, pos)) != std::string::npos) + { + ++count; + pos += pattern.size(); + } + return count; +} + +void init_unitcell(UnitCell& ucell) +{ + ucell.latName = "user_defined_lattice"; + ucell.lat0 = 10.0; + ucell.latvec.e11 = 1.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e33 = 1.0; + ucell.ntype = 1; + ucell.nat = 1; + ucell.atoms = new Atom[1]; + ucell.set_atom_flag = true; + ucell.atoms[0].label = "Si"; + ucell.atoms[0].na = 1; + ucell.atoms[0].nw = 2; + ucell.atoms[0].taud.resize(1); + ucell.atoms[0].taud[0] = ModuleBase::Vector3(0.0, 0.25, 0.5); +} + +void init_serial_orbitals(Parallel_Orbitals& pv) +{ + pv.atom_begin_row.resize(2); + pv.atom_begin_col.resize(2); + pv.atom_begin_row[0] = 0; + pv.atom_begin_row[1] = 2; + pv.atom_begin_col[0] = 0; + pv.atom_begin_col[1] = 2; + pv.nrow = 2; + pv.ncol = 2; + pv.set_serial(2, 2); +} + +void fill_matrix(hamilt::HContainer& matrix, Parallel_Orbitals& pv, double* values) +{ + hamilt::AtomPair pair(0, 0, 0, 0, 0, &pv, values); + matrix.insert_pair(pair); +} + +bool starts_with(const std::string& text, const std::string& prefix) +{ + return text.find(prefix) == 0; +} +} // namespace + +TEST(WriteHsRCompatibility, FileNameHelpersKeepCurrentContract) +{ + EXPECT_EQ(ModuleIO::hsr_gen_fname("hrs", 0, true, -1), "hrs1_nao.csr"); + EXPECT_EQ(ModuleIO::hsr_gen_fname("hrs", 1, true, 0), "hrs2_nao.csr"); + EXPECT_EQ(ModuleIO::hsr_gen_fname("hrs", 1, false, 3), "hrs2g4_nao.csr"); + EXPECT_EQ(ModuleIO::hsr_gen_fname("srs", 0, false, 0), "srs1g1_nao.csr"); + EXPECT_EQ(ModuleIO::hsr_gen_fname("srs", 0, false, -1), "srs1_nao.csr"); + + EXPECT_EQ(ModuleIO::dhr_gen_fname("dhrx", 0, true, -1), "dhrxrs1_nao.csr"); + EXPECT_EQ(ModuleIO::dhr_gen_fname("dhrx", 0, false, 0), "dhrxrs1g1_nao.csr"); + EXPECT_EQ(ModuleIO::dhr_gen_fname("dsry", 1, false, 2), "dsryrs2g3_nao.csr"); + + EXPECT_EQ(ModuleIO::dmr_gen_fname(1, 0, true, -1), "dmrs1_nao.csr"); + EXPECT_EQ(ModuleIO::dmr_gen_fname(1, 1, false, 2), "dmrs2g3_nao.csr"); + EXPECT_EQ(ModuleIO::dmr_gen_fname(2, 0, true, 5), "dmrs1_nao.npz"); +} + +TEST(WriteHsRCompatibility, HContainerCsrHeaderKeepsCurrentFormat) +{ + const std::string filename = "write_hs_r_header_h.csr"; + std::remove(filename.c_str()); + + UnitCell ucell; + init_unitcell(ucell); + Parallel_Orbitals pv; + init_serial_orbitals(pv); + hamilt::HContainer matrix(&pv); + double values[4] = {1.0, 0.0, 0.5, 2.0}; + fill_matrix(matrix, pv, values); + + ModuleIO::write_hcontainer_csr(filename, &ucell, 5, &matrix, 0, 0, 1, "H"); + + const std::string output = read_file(filename); + EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 1 ---\n")); + EXPECT_THAT(output, testing::HasSubstr(" # print H matrix in real space H(R)\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1 # number of spin directions\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1 # spin index\n")); + EXPECT_THAT(output, testing::HasSubstr(" 2 # number of localized basis\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1 # number of Bravais lattice vector R\n")); + EXPECT_THAT(output, testing::HasSubstr(" user_defined_lattice\n")); + EXPECT_THAT(output, testing::HasSubstr(" Si\n")); + EXPECT_THAT(output, testing::HasSubstr(" Direct\n")); + EXPECT_THAT(output, testing::HasSubstr(" 0 0.25 0.5\n")); + EXPECT_THAT(output, testing::HasSubstr(" # CSR Format #\n")); + EXPECT_THAT(output, testing::HasSubstr(" 0 0 0 3\n")); + EXPECT_THAT(output, testing::HasSubstr(" # CSR values\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1.00000e+00 5.00000e-01 2.00000e+00")); + + std::remove(filename.c_str()); +} + +TEST(WriteHsRCompatibility, HContainerCsrAppendKeepsCurrentStepSections) +{ + const std::string filename = "write_hs_r_append_s.csr"; + std::remove(filename.c_str()); + + UnitCell ucell; + init_unitcell(ucell); + Parallel_Orbitals pv; + init_serial_orbitals(pv); + hamilt::HContainer matrix(&pv); + double values[4] = {1.0, 0.0, 0.0, 1.0}; + fill_matrix(matrix, pv, values); + + ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 0, 0, 1, "S"); + ModuleIO::write_hcontainer_csr(filename, &ucell, 4, &matrix, 1, 0, 1, "S"); + + const std::string output = read_file(filename); + EXPECT_EQ(count_substr(output, " --- Ionic Step "), 2); + EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 1 ---\n")); + EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 2 ---\n")); + EXPECT_EQ(count_substr(output, " # print S matrix in real space S(R)\n"), 2); + + std::remove(filename.c_str()); +} + +TEST(WriteHsRCompatibility, DmrCsrHeaderKeepsCurrentFormat) +{ + std::string filename = "write_hs_r_header_dmr.csr"; + std::remove(filename.c_str()); + + UnitCell ucell; + init_unitcell(ucell); + Parallel_Orbitals pv; + init_serial_orbitals(pv); + hamilt::HContainer matrix(&pv); + double values[4] = {0.25, 0.0, 0.0, 0.75}; + fill_matrix(matrix, pv, values); + + ModuleIO::write_dmr_csr(filename, &ucell, 4, &matrix, 0, 1, 2); + + const std::string output = read_file(filename); + EXPECT_THAT(output, testing::HasSubstr(" --- Ionic Step 1 ---\n")); + EXPECT_THAT(output, testing::HasSubstr(" # print density matrix in real space DM(R)\n")); + EXPECT_THAT(output, testing::HasSubstr(" 2 # number of spin directions\n")); + EXPECT_THAT(output, testing::HasSubstr(" 2 # spin index\n")); + EXPECT_THAT(output, testing::HasSubstr(" 2 # number of localized basis\n")); + EXPECT_THAT(output, testing::HasSubstr(" 1 # number of Bravais lattice vector R\n")); + + std::remove(filename.c_str()); +} + +TEST(WriteHsRCompatibility, LegacySparseHeaderKeepsStepStyle) +{ + const std::string filename = "write_hs_r_legacy_s.csr"; + std::remove(filename.c_str()); + + GlobalV::DRANK = 0; + PARAM.sys.global_out_dir = "./"; + PARAM.sys.nlocal = 2; + + Parallel_Orbitals pv; + init_serial_orbitals(pv); + const Abfs::Vector3_Order r_vector(0, 0, 0); + std::set> all_R_coor; + all_R_coor.insert(r_vector); + std::map, std::map>> sparse_matrix; + sparse_matrix[r_vector][0][1] = 0.5; + sparse_matrix[r_vector][1][1] = 1.5; + + ModuleIO::save_sparse(sparse_matrix, all_R_coor, 1e-10, false, filename, pv, "S", 0, false); + + const std::string output = read_file(filename); + EXPECT_TRUE(starts_with(output, "STEP: 0\n")); + EXPECT_THAT(output, testing::HasSubstr("Matrix Dimension of S(R): 2\n")); + EXPECT_THAT(output, testing::HasSubstr("Matrix number of S(R): 1\n")); + EXPECT_THAT(output, testing::HasSubstr("0 0 0 2\n")); + + std::remove(filename.c_str()); +} + +TEST(WriteHsRCompatibility, LegacySparseBinaryHeaderWritesConcreteStep) +{ + const std::string filename = "write_hs_r_legacy_binary_s.csr"; + std::remove(filename.c_str()); + + GlobalV::DRANK = 0; + PARAM.sys.global_out_dir = "./"; + PARAM.sys.nlocal = 2; + + Parallel_Orbitals pv; + init_serial_orbitals(pv); + const Abfs::Vector3_Order r_vector(0, 0, 0); + std::set> all_R_coor; + all_R_coor.insert(r_vector); + std::map, std::map>> sparse_matrix; + sparse_matrix[r_vector][0][1] = 0.5; + sparse_matrix[r_vector][1][1] = 1.5; + + ModuleIO::save_sparse(sparse_matrix, all_R_coor, 1e-10, true, filename, pv, "S", 3, false); + + const std::vector header_and_r = read_binary_ints(filename, 7); + EXPECT_THAT(header_and_r, testing::ElementsAre(3, 2, 1, 0, 0, 0, 2)); + + std::remove(filename.c_str()); +} + +TEST(WriteHsRCompatibility, HeaderStyleSamplesRemainDistinct) +{ + EXPECT_TRUE(starts_with(" --- Ionic Step 1 ---", " --- Ionic Step")); + EXPECT_TRUE(starts_with("STEP: 0", "STEP:")); + EXPECT_TRUE(starts_with("IONIC_STEP: 1", "IONIC_STEP:")); +}