diff --git a/cmake/FindMKL.cmake b/cmake/FindMKL.cmake index 3680b8a92c..1036a90b87 100644 --- a/cmake/FindMKL.cmake +++ b/cmake/FindMKL.cmake @@ -12,6 +12,10 @@ find_path(MKL_INCLUDE mkl_service.h HINTS ${MKLROOT}/include) find_library(MKL_INTEL NAMES mkl_intel_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64) find_library(MKL_INTEL_THREAD NAMES mkl_intel_thread HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64) find_library(MKL_CORE NAMES mkl_core HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64) +find_library(MKL_IOMP5 NAMES iomp5 + HINTS ENV CMPLR_ROOT + PATH_SUFFIXES lib lib/intel64 linux/compiler/lib/intel64_lin +) if(ENABLE_MPI) find_library(MKL_SCALAPACK NAMES mkl_scalapack_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64) find_library(MKL_BLACS_INTELMPI NAMES mkl_blacs_intelmpi_lp64 HINTS ${MKLROOT}/lib ${MKLROOT}/lib/intel64) @@ -58,6 +62,11 @@ if(MKL_FOUND) IMPORTED_LOCATION "${MKL_BLACS_INTELMPI}" INTERFACE_INCLUDE_DIRECTORIES "${MKL_INCLUDE}") endif() + if(MKL_IOMP5 AND NOT TARGET MKL::IOMP5) + add_library(MKL::IOMP5 UNKNOWN IMPORTED) + set_target_properties(MKL::IOMP5 PROPERTIES + IMPORTED_LOCATION "${MKL_IOMP5}") + endif() add_library(MKL::MKL INTERFACE IMPORTED) if (ENABLE_MPI) set_property(TARGET MKL::MKL PROPERTY @@ -74,6 +83,10 @@ if(MKL_FOUND) "-Wl,--end-group" ) endif() + if(TARGET MKL::IOMP5) + set_property(TARGET MKL::MKL APPEND PROPERTY + INTERFACE_LINK_LIBRARIES MKL::IOMP5) + endif() endif() if(ENABLE_MPI) diff --git a/docs/advanced/elec_properties/density_matrix.md b/docs/advanced/elec_properties/density_matrix.md index bc293a617f..15b7ce9f44 100644 --- a/docs/advanced/elec_properties/density_matrix.md +++ b/docs/advanced/elec_properties/density_matrix.md @@ -38,3 +38,15 @@ The following line is dimension of the density matrix, and the rest lines are th The examples can be found in [examples/density_matrix](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/density_matrix) - Note: now this function is valid only for LCAO gamma only calcualtion. + +## Real-space Density Matrix (CSR format) + +ABACUS can also output the real-space density matrix DM(R) in CSR (Compressed Sparse Row) format by setting: +``` +out_dmr 1 +``` +After the calculation, the density matrix files are written to `OUT.${suffix}/`: +- `nspin=1`: `dmrs1_nao.csr` +- `nspin=2` (spin-polarized): `dmrs1_nao.csr` (spin-up) and `dmrs2_nao.csr` (spin-down) + +These files can be used to restart calculations by setting `init_chg dm` in the INPUT file together with `read_file_dir` pointing to the directory containing the CSR files. This is supported for both `nspin=1` and `nspin=2`. diff --git a/docs/advanced/scf/initialization.md b/docs/advanced/scf/initialization.md index d3ad5e48ed..8332c85b7f 100644 --- a/docs/advanced/scf/initialization.md +++ b/docs/advanced/scf/initialization.md @@ -1,15 +1,17 @@ # Initializing SCF Good initializing would abate the number of iteration steps in SCF. -Charge density should be initialed for constructing the initial hamiltonian operator. +Charge density should be initialed for constructing the initial hamiltonian operator. In PW basis, wavefunction should be initialized for iterate diagonalization method. In LCAO basis, wavefunction can be read to calculate initial charge density. The wavefunction itself does not have to be initialized. ## Charge Density `init_chg` is used for choosing the method of charge density initialization. - - `atomic` : initial charge density by atomic charge density from pseudopotential file under keyword `PP_RHOATOM` + - `atomic` : initial charge density by atomic charge density from pseudopotential file under keyword `PP_RHOATOM` - `file` : initial charge density from files produced by previous calculations with [`out_chg 1`](../elec_properties/charge.md). - `auto`: Abacus first attempts to read the density from a file; if not found, it defaults to using atomic density. + - `dm` (LCAO only): initial charge density from density matrix files in CSR format. For `nspin=1`, reads `dmrs1_nao.csr`. For `nspin=2` (spin-polarized), reads both `dmrs1_nao.csr` (spin-up) and `dmrs2_nao.csr` (spin-down). These files are generated by previous calculations with [`out_dmr 1`](../elec_properties/density_matrix.md). This method is particularly useful for restarting spin-polarized calculations. + - `hr` (LCAO only): initial charge density from Hamiltonian matrix files in CSR format. The Hamiltonian is read from file, then diagonalized to obtain wavefunctions and charge density. For `nspin=1`, reads `hrs1_nao.csr`. For `nspin=2` (spin-polarized), reads both `hrs1_nao.csr` (spin-up) and `hrs2_nao.csr` (spin-down). These files are generated by previous calculations with [`out_mat_hs2 1`](../input_files/input-main.md). ## Wave function `init_wfc` is used for choosing the method of wavefunction coefficient initialization. diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 0e7cbef3ec..f8cecf6805 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -170,31 +170,21 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->dmat.dm->init_DMR(*hamilt_lcao->getHR()); // 13.1) decide the strategy for initializing DMR and HR - if(istep == 0)//if the first scf step, readin DMR from file, + if(istep == 0)//if the first scf step, readin DMR from file, { //calculate or readin the density matrix DMR if(PARAM.inp.init_chg == "dm") { - //! 13.1.1) init density matrix from file - std::string dmfile = PARAM.globalv.global_readin_dir + "dmrs1_nao.csr"; - LCAO_domain::init_dm_from_file(dmfile, this->dmat, ucell, &(this->pv)); - GlobalV::ofs_running << " Read density matrix (real space) from file: " - << dmfile << std::endl; + //! 13.1.1) init charge density from density matrix file + LCAO_domain::init_chg_dm(PARAM.globalv.global_readin_dir, PARAM.inp.nspin, + this->dmat, ucell, &(this->pv), this->pelec->charge); } if(PARAM.inp.init_chg == "hr") { - //! 13.1.2) init HR from file - std::string hrfile = PARAM.globalv.global_readin_dir + "hrs1_nao.csr"; - LCAO_domain::init_hr_from_file( - hrfile, - dynamic_cast*>(this->p_hamilt)->getHR(), - ucell, - &(this->pv) - ); - this->p_hamilt->refresh(false); - hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); - hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, - this->chr, PARAM.inp.nspin, 0); + //! 13.1.2) init charge density from Hamiltonian matrix file + LCAO_domain::init_chg_hr(PARAM.globalv.global_readin_dir, PARAM.inp.nspin, + this->p_hamilt, ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm, + this->chr, PARAM.inp.ks_solver); } } else //if not, use the DMR calculated from last step @@ -204,11 +194,6 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros this->dmat.dm->cal_DMR(); } - // 13.2 if init_chg = "dm", then calculate rho from readin DMR before init_scf - if(PARAM.inp.init_chg == "dm") - { - LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, this->pelec->charge, true); - } // 13.2) init_scf, should be before_scf? mohan add 2025-03-10 this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt index 2bdbfc0b45..844da5dc84 100644 --- a/source/source_lcao/CMakeLists.txt +++ b/source/source_lcao/CMakeLists.txt @@ -64,6 +64,10 @@ if(ENABLE_LCAO) add_coverage(hamilt_lcao) endif() + IF (BUILD_TESTING) + add_subdirectory(test) + endif() + endif() diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp index 9d30b073a0..2b38f94f93 100644 --- a/source/source_lcao/LCAO_set.cpp +++ b/source/source_lcao/LCAO_set.cpp @@ -4,6 +4,9 @@ #include "source_io/module_wf/read_wfc_nao.h" // use read_wfc_nao #include "source_estate/elecstate_tools.h" // use fixed_weights #include "source_lcao/module_hcontainer/read_hcontainer.h" +#include "source_lcao/rho_tau_lcao.h" // use dm2rho +#include "source_lcao/hamilt_lcao.h" // use HamiltLCAO for init_chg_hr +#include "source_hsolver/hsolver_lcao.h" // use HSolverLCAO for init_chg_hr template void LCAO_domain::set_psi_occ_dm_chg( @@ -94,20 +97,46 @@ void LCAO_domain::set_pot( template void LCAO_domain::init_dm_from_file( - const std::string dmfile, + const std::string& readin_dir, + const int nspin, LCAO_domain::Setup_DM& dmat, const UnitCell& ucell, const Parallel_Orbitals* pv) { ModuleBase::TITLE("LCAO_domain", "init_dm_from_file"); - hamilt::HContainer* dm_container = dmat.dm->get_DMR_vector()[0]; - hamilt::Read_HContainer reader_dm( - dm_container, - dmfile, - PARAM.globalv.nlocal, - &ucell - ); - reader_dm.read(); + const int nspin_dm = (nspin == 2) ? 2 : 1; + for (int is = 0; is < nspin_dm; ++is) + { + const std::string dmfile = readin_dir + "/dmrs" + std::to_string(is + 1) + "_nao.csr"; + hamilt::HContainer* dm_container = dmat.dm->get_DMR_vector()[is]; + hamilt::Read_HContainer reader_dm( + dm_container, + dmfile, + PARAM.globalv.nlocal, + &ucell + ); + reader_dm.read(); + } + return; +} + +template +void LCAO_domain::init_chg_dm( + const std::string& readin_dir, + const int nspin, + LCAO_domain::Setup_DM& dmat, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + Charge* chr) +{ + ModuleBase::TITLE("LCAO_domain", "init_chg_dm"); + + // Step 1: Read density matrix from file + LCAO_domain::init_dm_from_file(readin_dir, nspin, dmat, ucell, pv); + + // Step 2: Convert density matrix to charge density + LCAO_domain::dm2rho(dmat.dm->get_DMR_vector(), nspin, chr, true); + return; } @@ -130,6 +159,57 @@ void LCAO_domain::init_hr_from_file( return; } +template +void LCAO_domain::init_chg_hr( + const std::string& readin_dir, + const int nspin, + hamilt::Hamilt* p_hamilt, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + psi::Psi& psi, + elecstate::ElecState* pelec, + elecstate::DensityMatrix& dm, + Charge& chr, + const std::string& ks_solver) +{ + ModuleBase::TITLE("LCAO_domain", "init_chg_hr"); + + auto* hamilt_lcao = dynamic_cast*>(p_hamilt); + if (!hamilt_lcao) + { + ModuleBase::WARNING_QUIT("LCAO_domain::init_chg_hr", "p_hamilt is not HamiltLCAO"); + } + + // Step 1: Read HR from file(s) + if (nspin == 2) + { + // nspin=2: load spin-up into first half of hRS2, spin-down into second half + const std::string hrfile_up = readin_dir + "/hrs1_nao.csr"; + LCAO_domain::init_hr_from_file(hrfile_up, hamilt_lcao->getHR(), ucell, pv); + + // switch hR data pointer to spin-down half, then read hrs2 + auto& hRS2 = hamilt_lcao->getHRS2(); + hamilt_lcao->getHR()->allocate(hRS2.data() + hRS2.size() / 2, 0); + const std::string hrfile_down = readin_dir + "/hrs2_nao.csr"; + LCAO_domain::init_hr_from_file(hrfile_down, hamilt_lcao->getHR(), ucell, pv); + + // restore hR to spin-up half + hamilt_lcao->getHR()->allocate(hRS2.data(), 0); + } + else + { + const std::string hrfile = readin_dir + "/hrs1_nao.csr"; + LCAO_domain::init_hr_from_file(hrfile, hamilt_lcao->getHR(), ucell, pv); + } + + // Step 2: Mark HR as loaded from file (skip operator recalculation) + p_hamilt->refresh(false); + + // Step 3: Diagonalize to get wavefunctions and charge density + hsolver::HSolverLCAO hsolver_lcao_obj(pv, ks_solver); + hsolver_lcao_obj.solve(p_hamilt, psi, pelec, dm, chr, nspin, 0); +} + template void LCAO_domain::set_psi_occ_dm_chg( @@ -183,16 +263,33 @@ template void LCAO_domain::set_pot>( const Input_para &inp); template void LCAO_domain::init_dm_from_file( - const std::string dmfile, + const std::string& readin_dir, + const int nspin, LCAO_domain::Setup_DM& dmat, const UnitCell& ucell, const Parallel_Orbitals* pv); template void LCAO_domain::init_dm_from_file>( - const std::string dmfile, + const std::string& readin_dir, + const int nspin, LCAO_domain::Setup_DM>& dmat, const UnitCell& ucell, const Parallel_Orbitals* pv); +template void LCAO_domain::init_chg_dm( + const std::string& readin_dir, + const int nspin, + LCAO_domain::Setup_DM& dmat, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + Charge* chr); +template void LCAO_domain::init_chg_dm>( + const std::string& readin_dir, + const int nspin, + LCAO_domain::Setup_DM>& dmat, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + Charge* chr); + template void LCAO_domain::init_hr_from_file( const std::string hrfile, hamilt::HContainer* hmat, @@ -203,3 +300,37 @@ template void LCAO_domain::init_hr_from_file>( hamilt::HContainer>* hmat, const UnitCell& ucell, const Parallel_Orbitals* pv); + +template void LCAO_domain::init_chg_hr( + const std::string& readin_dir, + const int nspin, + hamilt::Hamilt* p_hamilt, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + psi::Psi& psi, + elecstate::ElecState* pelec, + elecstate::DensityMatrix& dm, + Charge& chr, + const std::string& ks_solver); +template void LCAO_domain::init_chg_hr, double>( + const std::string& readin_dir, + const int nspin, + hamilt::Hamilt>* p_hamilt, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + psi::Psi>& psi, + elecstate::ElecState* pelec, + elecstate::DensityMatrix, double>& dm, + Charge& chr, + const std::string& ks_solver); +template void LCAO_domain::init_chg_hr, std::complex>( + const std::string& readin_dir, + const int nspin, + hamilt::Hamilt>* p_hamilt, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + psi::Psi>& psi, + elecstate::ElecState* pelec, + elecstate::DensityMatrix, double>& dm, + Charge& chr, + const std::string& ks_solver); diff --git a/source/source_lcao/LCAO_set.h b/source/source_lcao/LCAO_set.h index fc7d807a37..f086bc5027 100644 --- a/source/source_lcao/LCAO_set.h +++ b/source/source_lcao/LCAO_set.h @@ -5,6 +5,8 @@ #include "source_cell/klist.h" #include "source_psi/psi.h" #include "source_estate/elecstate.h" +#include "source_estate/module_dm/density_matrix.h" +#include "source_hamilt/hamilt.h" #include "source_lcao/setup_dm.h" #include "source_pw/module_pwdft/structure_factor.h" #include "source_basis/module_pw/pw_basis.h" @@ -55,14 +57,36 @@ void set_pot( /** * @brief read in DMR from file, and save it into dmat + * @param readin_dir directory containing dmrs*_nao.csr files + * @param nspin number of spin components (1 or 2) */ template void init_dm_from_file( - const std::string dmfile, + const std::string& readin_dir, + const int nspin, LCAO_domain::Setup_DM& dmat, const UnitCell& ucell, const Parallel_Orbitals* pv); +/** + * @brief initialize charge density from density matrix file (init_chg=dm) + * This function reads DMR from file and converts it to charge density + * @param readin_dir directory containing dmrs*_nao.csr files + * @param nspin number of spin components (1 or 2) + * @param dmat density matrix object + * @param ucell unit cell + * @param pv parallel orbitals + * @param chr charge density object + */ +template +void init_chg_dm( + const std::string& readin_dir, + const int nspin, + LCAO_domain::Setup_DM& dmat, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + Charge* chr); + /** * @brief read in HR from file, and save it into hmat */ @@ -72,6 +96,37 @@ void init_hr_from_file( hamilt::HContainer* hmat, const UnitCell& ucell, const Parallel_Orbitals* pv); + +/** + * @brief initialize charge density from Hamiltonian matrix file (init_chg=hr) + * Reads HR from file(s), diagonalizes to get wavefunctions, then computes charge density. + * For nspin=2, reads both hrs1_nao.csr (spin-up) and hrs2_nao.csr (spin-down) + * into the two halves of HamiltLCAO::hRS2. + * @tparam TK k-space type (double or complex) + * @tparam TR real-space type (double) + * @param readin_dir directory containing hrs*_nao.csr files + * @param nspin number of spin components + * @param p_hamilt pointer to Hamilt base class (will be dynamic_cast to HamiltLCAO) + * @param ucell unit cell + * @param pv parallel orbitals + * @param psi wave function object + * @param pelec electronic state + * @param dm density matrix + * @param chr charge density + * @param ks_solver solver method name + */ +template +void init_chg_hr( + const std::string& readin_dir, + const int nspin, + hamilt::Hamilt* p_hamilt, + const UnitCell& ucell, + const Parallel_Orbitals* pv, + psi::Psi& psi, + elecstate::ElecState* pelec, + elecstate::DensityMatrix& dm, + Charge& chr, + const std::string& ks_solver); } // end namespace #endif diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index 1b3d7195a6..3850f636bd 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -537,8 +537,10 @@ void HamiltLCAO::refresh(bool yes) this->refresh_times = 0; if (PARAM.inp.nspin == 2) { - ModuleBase::WARNING_QUIT("HamiltLCAO::refresh", - "When turning off the refresh flag, the nspin==2 case is not supported yet."); + // HR has been loaded from file into both halves of hRS2. + // Reset to spin-up; updateHk will switch pointers as needed. + this->current_spin = 0; + this->hR->allocate(this->hRS2.data(), 0); } } } diff --git a/source/source_lcao/hamilt_lcao.h b/source/source_lcao/hamilt_lcao.h index deca3c7c73..49b309b8c3 100644 --- a/source/source_lcao/hamilt_lcao.h +++ b/source/source_lcao/hamilt_lcao.h @@ -123,6 +123,9 @@ class HamiltLCAO : public Hamilt } #endif + /// get hRS2 buffer for NSPIN=2 case (spin-up in first half, spin-down in second half) + std::vector& getHRS2() { return this->hRS2; } + /// refresh the status of HR void refresh(bool yes) override; diff --git a/source/source_lcao/test/CMakeLists.txt b/source/source_lcao/test/CMakeLists.txt new file mode 100644 index 0000000000..cd2a887cc3 --- /dev/null +++ b/source/source_lcao/test/CMakeLists.txt @@ -0,0 +1,25 @@ +remove_definitions(-D__MLALGO) +remove_definitions(-D__CUDA) +remove_definitions(-D__ROCM) + +if(ENABLE_LCAO) +AddTest( + TARGET init_dm_from_file_test + LIBS parameter ${math_libs} base device + SOURCES test_init_dm_from_file.cpp tmp_mocks.cpp + ${ABACUS_SOURCE_DIR}/source_estate/module_dm/density_matrix.cpp + ${ABACUS_SOURCE_DIR}/source_estate/module_dm/density_matrix_io.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/base_matrix.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/hcontainer.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/atom_pair.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/read_hcontainer.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/func_transfer.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/func_folding.cpp + ${ABACUS_SOURCE_DIR}/source_lcao/module_hcontainer/transfer.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp + ${ABACUS_SOURCE_DIR}/source_io/module_output/sparse_matrix.cpp + ${ABACUS_SOURCE_DIR}/source_io/module_output/csr_reader.cpp + ${ABACUS_SOURCE_DIR}/source_io/module_output/file_reader.cpp + ${ABACUS_SOURCE_DIR}/source_io/module_output/output.cpp +) +endif() diff --git a/source/source_lcao/test/test_init_dm_from_file.cpp b/source/source_lcao/test/test_init_dm_from_file.cpp new file mode 100644 index 0000000000..ea23f82dc7 --- /dev/null +++ b/source/source_lcao/test/test_init_dm_from_file.cpp @@ -0,0 +1,252 @@ +#include +#include + +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#define private public +#include "source_estate/module_dm/density_matrix.h" +#include "source_lcao/module_hcontainer/hcontainer.h" +#include "source_lcao/module_hcontainer/read_hcontainer.h" +#include "source_lcao/setup_dm.h" +#include "source_cell/klist.h" +#undef private + +/************************************************ + * unit test of init_dm_from_file (nspin=1 & nspin=2) + ***********************************************/ + +// Small test system: 2 atoms, 13 orbitals each => nlocal=26 +int test_size = 2; +int test_nw = 13; + +class InitDMFileTest : public testing::Test +{ + protected: + Parallel_Orbitals* paraV; + UnitCell ucell; + int nlocal; + + void SetUp() override + { + nlocal = test_size * test_nw; + + // set up a unitcell + ucell.ntype = 1; + ucell.nat = test_size; + ucell.atoms = new Atom[ucell.ntype]; + ucell.iat2it = new int[ucell.nat]; + ucell.iat2ia = new int[ucell.nat]; + ucell.atoms[0].tau.resize(ucell.nat); + ucell.itia2iat.create(ucell.ntype, ucell.nat); + for (int iat = 0; iat < ucell.nat; iat++) + { + ucell.iat2it[iat] = 0; + ucell.iat2ia[iat] = iat; + ucell.atoms[0].tau[iat] = ModuleBase::Vector3(0.0, 0.0, 0.0); + ucell.itia2iat(0, iat) = iat; + } + ucell.atoms[0].na = test_size; + ucell.atoms[0].nw = test_nw; + ucell.atoms[0].iw2l.resize(test_nw, 0); + ucell.atoms[0].iw2m.resize(test_nw, 0); + ucell.atoms[0].iw2n.resize(test_nw, 0); + ucell.set_iat2iwt(1); + + // set up parallel orbitals (serial mode) + paraV = new Parallel_Orbitals(); + paraV->set_serial(nlocal, nlocal); + paraV->set_atomic_trace(ucell.get_iat2iwt(), ucell.nat, nlocal); + } + + void TearDown() override + { + delete paraV; + delete[] ucell.atoms; + } + + /// Write a minimal CSR file with a diagonal matrix at R=(0,0,0) + void write_test_csr(const std::string& filename, double scale) + { + std::ofstream ofs(filename); + ofs << "IONIC_STEP: 1" << std::endl; + ofs << "Matrix Dimension of DM(R): " << nlocal << std::endl; + ofs << "Matrix number of DM(R): 1" << std::endl; + + // R coordinate header: rx ry rz nnz + int nnz = nlocal; // diagonal + ofs << "0 0 0 " << nnz << std::endl; + + // values line + for (int i = 0; i < nlocal; i++) + { + if (i > 0) ofs << " "; + ofs << std::scientific << std::setprecision(8) << scale * (i + 1) * 0.01; + } + ofs << std::endl; + + // column indices line + for (int i = 0; i < nlocal; i++) + { + if (i > 0) ofs << " "; + ofs << i; + } + ofs << std::endl; + + // row pointers line (CSR format: each row has exactly 1 element on diagonal) + for (int i = 0; i <= nlocal; i++) + { + if (i > 0) ofs << " "; + ofs << i; + } + ofs << std::endl; + + ofs.close(); + } + + /// Create DensityMatrix with given nspin and initialize DMR from an HContainer template + elecstate::DensityMatrix* create_dm(int nspin) + { + K_Vectors kv; + int nks = (nspin == 2) ? 2 : 1; // gamma_only: nk=1 per spin + kv.set_nks(nks * (nspin == 2 ? 2 : 1)); + kv.kvec_d.resize(kv.get_nks()); + + int nspin_dm = (nspin == 2) ? 2 : 1; + auto* dm = new elecstate::DensityMatrix( + paraV, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); + + // Create a template HContainer and init DMR from it + hamilt::HContainer tmp_HR(paraV); + // Add atom pairs for all atom-atom combinations at R=(0,0,0) + for (int i = 0; i < ucell.nat; i++) + { + for (int j = 0; j < ucell.nat; j++) + { + hamilt::AtomPair ap(i, j, 0, 0, 0, paraV); + tmp_HR.insert_pair(ap); + } + } + tmp_HR.allocate(nullptr, true); + dm->init_DMR(tmp_HR); + return dm; + } +}; + +TEST_F(InitDMFileTest, Nspin1_ReadSingleFile) +{ + // Create test directory and CSR file + system("mkdir -p ./test_dm_dir"); + write_test_csr("./test_dm_dir/dmrs1_nao.csr", 1.0); + + // Create DM with nspin=1 + auto* dm = create_dm(1); + ASSERT_EQ(dm->_DMR.size(), 1); + + // Read from file using Read_HContainer (same as init_dm_from_file does) + hamilt::HContainer* dmr0 = dm->get_DMR_vector()[0]; + hamilt::Read_HContainer reader(dmr0, "./test_dm_dir/dmrs1_nao.csr", nlocal, &ucell); + reader.read(); + + // Verify DMR[0] has data + EXPECT_GT(dmr0->size_atom_pairs(), 0); + + // Check diagonal element (0,0) at R=(0,0,0) is non-zero + auto* ap = dmr0->find_pair(0, 0); + ASSERT_NE(ap, nullptr); + bool has_nonzero = false; + for (int i = 0; i < ap->get_size(); i++) + { + if (std::abs(ap->get_pointer()[i]) > 1e-15) + { + has_nonzero = true; + break; + } + } + EXPECT_TRUE(has_nonzero); + + delete dm; + system("rm -rf ./test_dm_dir"); +} + +TEST_F(InitDMFileTest, Nspin2_ReadTwoFiles) +{ + // Create test directory and two CSR files with different scale factors + system("mkdir -p ./test_dm_dir"); + write_test_csr("./test_dm_dir/dmrs1_nao.csr", 1.0); // spin-up + write_test_csr("./test_dm_dir/dmrs2_nao.csr", 0.5); // spin-down + + // Create DM with nspin=2 + auto* dm = create_dm(2); + ASSERT_EQ(dm->_DMR.size(), 2); + + // Read spin-up + hamilt::HContainer* dmr0 = dm->get_DMR_vector()[0]; + hamilt::Read_HContainer reader0(dmr0, "./test_dm_dir/dmrs1_nao.csr", nlocal, &ucell); + reader0.read(); + + // Read spin-down + hamilt::HContainer* dmr1 = dm->get_DMR_vector()[1]; + hamilt::Read_HContainer reader1(dmr1, "./test_dm_dir/dmrs2_nao.csr", nlocal, &ucell); + reader1.read(); + + // Verify both DMR components have data + EXPECT_GT(dmr0->size_atom_pairs(), 0); + EXPECT_GT(dmr1->size_atom_pairs(), 0); + + // Verify spin-up and spin-down have different values (scale 1.0 vs 0.5) + auto* ap0 = dmr0->find_pair(0, 0); + auto* ap1 = dmr1->find_pair(0, 0); + ASSERT_NE(ap0, nullptr); + ASSERT_NE(ap1, nullptr); + + bool values_differ = false; + int check_size = std::min(ap0->get_size(), ap1->get_size()); + for (int i = 0; i < check_size; i++) + { + double v0 = ap0->get_pointer()[i]; + double v1 = ap1->get_pointer()[i]; + if (std::abs(v0) > 1e-15 && std::abs(v0 - v1) > 1e-15) + { + values_differ = true; + // spin-down should be ~half of spin-up + EXPECT_NEAR(v1 / v0, 0.5, 1e-6); + break; + } + } + EXPECT_TRUE(values_differ); + + delete dm; + system("rm -rf ./test_dm_dir"); +} + +TEST_F(InitDMFileTest, Nspin2_DMRVectorSize) +{ + // Verify that nspin=2 creates exactly 2 DMR components + auto* dm = create_dm(2); + EXPECT_EQ(dm->_DMR.size(), 2); + EXPECT_NE(dm->_DMR[0], nullptr); + EXPECT_NE(dm->_DMR[1], nullptr); + delete dm; +} + +TEST_F(InitDMFileTest, Nspin1_DMRVectorSize) +{ + // Verify that nspin=1 creates exactly 1 DMR component + auto* dm = create_dm(1); + EXPECT_EQ(dm->_DMR.size(), 1); + EXPECT_NE(dm->_DMR[0], nullptr); + delete dm; +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} diff --git a/source/source_lcao/test/tmp_mocks.cpp b/source/source_lcao/test/tmp_mocks.cpp new file mode 100644 index 0000000000..515d80538b --- /dev/null +++ b/source/source_lcao/test/tmp_mocks.cpp @@ -0,0 +1,58 @@ + +#include "source_cell/unitcell.h" +#include "source_cell/module_neighbor/sltk_grid_driver.h" + +// constructor of Atom +Atom::Atom() {} +Atom::~Atom() {} + +Atom_pseudo::Atom_pseudo() {} +Atom_pseudo::~Atom_pseudo() {} + +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} +LCAO_Orbitals::LCAO_Orbitals() {} +LCAO_Orbitals::~LCAO_Orbitals() {} +#endif + +pseudo::pseudo() {} +pseudo::~pseudo() {} +SepPot::SepPot(){} +SepPot::~SepPot(){} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} + +// constructor of UnitCell +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} + +void UnitCell::set_iat2iwt(const int& npol_in) +{ + this->iat2iwt.resize(this->nat); + this->npol = npol_in; + int iat = 0; + int iwt = 0; + for (int it = 0; it < this->ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { + this->iat2iwt[iat] = iwt; + iwt += atoms[it].nw * this->npol; + ++iat; + } + } + return; +} + +// stub for Grid_Driver::Find_atom (used by density_matrix_io.cpp but not exercised in test) +void Grid_Driver::Find_atom(const UnitCell& ucell, + const ModuleBase::Vector3& tau, + const int& T, + const int& I, + AdjacentAtomInfo* adjs) const +{ +}