From 26072ee750c7f571a07880d4124662e87d66c976 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 2 Apr 2025 10:51:15 +0800 Subject: [PATCH 1/5] add input parameters for ldos --- docs/advanced/input_files/input-main.md | 17 ++++++++++++++++- source/module_io/read_input_item_output.cpp | 6 ++++++ .../module_io/read_input_item_postprocess.cpp | 12 ++++++++++++ source/module_io/test/read_input_ptest.cpp | 2 ++ source/module_io/test/support/INPUT | 2 ++ .../test_serial/read_input_item_test.cpp | 10 ++++++++++ source/module_parameter/input_parameter.h | 2 ++ 7 files changed, 50 insertions(+), 1 deletion(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index f8eec32d73d..81d873860b2 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -140,6 +140,7 @@ - [out\_wfc\_r](#out_wfc_r) - [out\_wfc\_lcao](#out_wfc_lcao) - [out\_dos](#out_dos) + - [out\_ldos](#out_ldos) - [out\_band](#out_band) - [out\_proj\_band](#out_proj_band) - [out\_stru](#out_stru) @@ -174,6 +175,7 @@ - [dos\_emin\_ev](#dos_emin_ev) - [dos\_emax\_ev](#dos_emax_ev) - [dos\_nche](#dos_nche) + - [stm\_bias](#stm_bias) - [NAOs](#naos) - [bessel\_nao\_ecut](#bessel_nao_ecut) - [bessel\_nao\_tolerence](#bessel_nao_tolerence) @@ -1706,6 +1708,12 @@ These variables are used to control the output of properties. - lcao-only: output the density of states (DOS) and the projected density of states (PDOS) - **Default**: 0 +### out_ldos + +- **Type**: Boolean +- **Description**: Whether to output the local density of states for given bias in cube file format, which is controlled by [stm_bias](#stm_bias). +- **Default**: False + ### out_band - **Type**: Boolean \[Integer\](optional) @@ -1955,9 +1963,16 @@ These variables are used to control the calculation of DOS. [Detailed introducti ### dos_nche - **Type**: Integer -The order of Chebyshev expansions when using Stochastic Density Functional Theory (SDFT) to calculate DOS. +- **Description**: The order of Chebyshev expansions when using Stochastic Density Functional Theory (SDFT) to calculate DOS. - **Default**: 100 +### stm_bias + +- **Type**: Real +- **Description**: The bias voltage used to calculate local density of states to simulate scanning tunneling microscope, see details in [out_ldos](#out_ldos). +- **Default**: 1.0 +- **Unit**: V + [back to top](#full-list-of-input-keywords) ## NAOs diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index 0614ca1d08e..abd923b2b29 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -143,6 +143,12 @@ void ReadInput::item_output() }; this->add_item(item); } + { + Input_Item item("out_ldos"); + item.annotation = "output local density of states"; + read_sync_bool(input.out_ldos); + this->add_item(item); + } { Input_Item item("out_mul"); item.annotation = "mulliken charge or not"; diff --git a/source/module_io/read_input_item_postprocess.cpp b/source/module_io/read_input_item_postprocess.cpp index 46c1e18b4f1..fec3cbf7fe6 100644 --- a/source/module_io/read_input_item_postprocess.cpp +++ b/source/module_io/read_input_item_postprocess.cpp @@ -53,6 +53,18 @@ void ReadInput::item_postprocess() read_sync_int(input.dos_nche); this->add_item(item); } + { + Input_Item item("stm_bias"); + item.annotation = "bias voltage used to calculate ldos"; + read_sync_double(input.stm_bias); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.out_ldos && para.input.stm_bias == 0.0) + { + ModuleBase::WARNING_QUIT("ReadInput", "a nonzero stm_bias is required for ldos calculation"); + } + }; + this->add_item(item); + } // Electronic Conductivity { diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index d551f4f5fc5..f72796268c1 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -201,6 +201,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.out_wfc_pw, 0); EXPECT_EQ(param.inp.out_wfc_r, 0); EXPECT_EQ(param.inp.out_dos, 0); + EXPECT_EQ(param.inp.out_ldos, true); EXPECT_EQ(param.inp.out_band[0], 0); EXPECT_EQ(param.inp.out_band[1], 8); EXPECT_EQ(param.inp.out_proj_band, 0); @@ -220,6 +221,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_DOUBLE_EQ(param.inp.dos_edelta_ev, 0.01); EXPECT_DOUBLE_EQ(param.inp.dos_scale, 0.01); EXPECT_DOUBLE_EQ(param.inp.dos_sigma, 0.07); + EXPECT_DOUBLE_EQ(param.inp.stm_bias, 2.0); EXPECT_FALSE(param.inp.out_element_info); EXPECT_DOUBLE_EQ(param.inp.lcao_ecut, 20); EXPECT_DOUBLE_EQ(param.inp.lcao_dk, 0.01); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 944e46ca2f8..f168ee2eab9 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -65,6 +65,7 @@ out_pot 2 #output realspace potential out_wfc_pw 0 #output wave functions out_wfc_r 0 #output wave functions in realspace out_dos 0 #output energy and dos +out_ldos True #output local density of states out_band 0 #output energy and band structure out_proj_band FaLse #output projected band structure restart_save f #print to disk every step for restart @@ -179,6 +180,7 @@ dos_edelta_ev 0.01 #delta energy for dos dos_scale 0.01 #scale dos range by dos_sigma 0.07 #gauss b coefficeinet(default=0.07) dos_nche 100 #orders of Chebyshev expansions for dos +stm_bias 2.0 #bias voltage used to calculate ldos #Parameters (9.Molecular dynamics) md_type nvt #choose ensemble diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp index d77cd1a3e96..21d3e3f3c25 100644 --- a/source/module_io/test_serial/read_input_item_test.cpp +++ b/source/module_io/test_serial/read_input_item_test.cpp @@ -391,6 +391,16 @@ TEST_F(InputTest, Item_test) output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("NOTICE")); } + { // stm_bias + auto it = find_label("stm_bias", readinput.input_lists); + param.input.out_ldos = true; + param.input.stm_bias = 0.0; + + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + } { // out_band auto it = find_label("out_band", readinput.input_lists); it->second.str_values = {"1"}; diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 77ee42a2ed2..8486ef255ba 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -363,6 +363,7 @@ struct Input_para int printe = 0; ///< Print out energy for each band for every printe step, default is scf_nmax std::vector out_band = {0, 8}; ///< band calculation pengfei 2014-10-13 int out_dos = 0; ///< dos calculation. mohan add 20090909 + bool out_ldos = false; ///< ldos calculation bool out_mul = false; ///< qifeng add 2019-9-10 bool out_proj_band = false; ///< projected band structure calculation jiyy add 2022-05-11 std::string out_level = "ie"; ///< control the output information. @@ -412,6 +413,7 @@ struct Input_para double dos_scale = 0.01; double dos_sigma = 0.07; ///< pengfei 2014-10-13 int dos_nche = 100; ///< orders of Chebyshev expansions for dos + double stm_bias = 1.0; ///< bias voltage for STM bool cal_cond = false; ///< calculate electronic conductivities double cond_che_thr = 1e-8; ///< control the error of Chebyshev expansions From 30219e5e2caa4d1d60723e3793ce4025cf25dd7d Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 2 Apr 2025 16:48:41 +0800 Subject: [PATCH 2/5] Feature: ldos in pw basis set --- source/Makefile.Objects | 1 + source/module_elecstate/CMakeLists.txt | 1 + source/module_elecstate/cal_ldos.cpp | 40 +++++++++++++++++++ source/module_elecstate/elecstate.h | 9 +++++ source/module_elecstate/elecstate_pw.h | 3 ++ source/module_elecstate/fp_energy.cpp | 12 +++--- source/module_elecstate/fp_energy.h | 2 +- .../test/elecstate_fp_energy_test.cpp | 17 ++++---- source/module_esolver/esolver_ks_pw.cpp | 6 +++ source/module_esolver/esolver_of.cpp | 8 ++-- source/module_io/out_ldos.h | 27 +++++++++++++ 11 files changed, 105 insertions(+), 21 deletions(-) create mode 100644 source/module_elecstate/cal_ldos.cpp create mode 100644 source/module_io/out_ldos.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9b690c5fdf8..f99cf3613fd 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -235,6 +235,7 @@ OBJS_ELECSTAT=elecstate.o\ H_TDDFT_pw.o\ pot_xc.o\ cal_ux.o\ + cal_ldos.o\ read_orb.o\ cal_nelec_nband.o\ read_pseudo.o\ diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 0bdd120879e..15e474789de 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -32,6 +32,7 @@ list(APPEND objects magnetism.cpp occupy.cpp cal_ux.cpp + cal_ldos.cpp read_orb.cpp cal_nelec_nband.cpp read_pseudo.cpp diff --git a/source/module_elecstate/cal_ldos.cpp b/source/module_elecstate/cal_ldos.cpp new file mode 100644 index 00000000000..8ef35a6f954 --- /dev/null +++ b/source/module_elecstate/cal_ldos.cpp @@ -0,0 +1,40 @@ +#include "elecstate_pw.h" + +namespace elecstate +{ +template +void ElecStatePW::cal_ldos(const psi::Psi& psi, std::vector& ldos) const +{ + // energy range for ldos (efermi as reference) + const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0; + const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0; + + ldos.resize(this->basis->nrxx); + std::vector wfcr(this->basis->nrxx); + + for (int ik = 0; ik < this->klist->get_nks(); ++ik) + { + psi.fix_k(ik); + double efermi = this->eferm.get_efval(this->klist->isk[ik]); + int nbands = psi.get_nbands(); + + for (int ib = 0; ib < nbands; ib++) + { + this->basis->recip2real(&psi(ib, 0), wfcr.data(), ik); + double eigenval = (ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV; + if (eigenval >= emin && eigenval <= emax) + { + for (int ir = 0; ir < basis->nrxx; ir++) + ldos[ir] += this->klist->wk[ik] * norm(wfcr[ir]); + } + } + } +} + +template class ElecStatePW, base_device::DEVICE_CPU>; +template class ElecStatePW, base_device::DEVICE_CPU>; +#if ((defined __CUDA) || (defined __ROCM)) +template class ElecStatePW, base_device::DEVICE_GPU>; +template class ElecStatePW, base_device::DEVICE_GPU>; +#endif +} // namespace elecstate diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index 5757ab696e5..e83e714f381 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -61,6 +61,15 @@ class ElecState return; } + virtual void cal_ldos(const psi::Psi>& psi, std::vector& ldos) const + { + return; + } + virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const + { + return; + } + // update charge density for next scf step // in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge // density rho among these rho, diff --git a/source/module_elecstate/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h index 679b9b712cb..da12f05a5c7 100644 --- a/source/module_elecstate/elecstate_pw.h +++ b/source/module_elecstate/elecstate_pw.h @@ -35,6 +35,9 @@ class ElecStatePW : public ElecState virtual void cal_tau(const psi::Psi& psi); + //! calculate local density of states + virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const; + //! calculate becsum for uspp void cal_becsum(const psi::Psi& psi); diff --git a/source/module_elecstate/fp_energy.cpp b/source/module_elecstate/fp_energy.cpp index 6855ef17c40..4920d66e2a4 100644 --- a/source/module_elecstate/fp_energy.cpp +++ b/source/module_elecstate/fp_energy.cpp @@ -96,22 +96,22 @@ void fenergy::print_all() const std::cout << " total= " << etot << std::endl; } -/// @brief get the reference of fermi of a specific spin +/// @brief set efermi of a specific spin /// @param is SPIN -/// @return a reference of fermi(is) -double& efermi::get_ef(const int& is) +/// @param ef_in fermi(is) +void efermi::set_efval(const int& is, const double& ef_in) { if (!two_efermi) { - return this->ef; + this->ef = ef_in; } else if (is == 0) { - return this->ef_up; + this->ef_up = ef_in; } else if (is == 1) { - return this->ef_dw; + this->ef_dw = ef_in; } else { diff --git a/source/module_elecstate/fp_energy.h b/source/module_elecstate/fp_energy.h index 17fad4c8beb..9d1e9c9fed7 100644 --- a/source/module_elecstate/fp_energy.h +++ b/source/module_elecstate/fp_energy.h @@ -64,7 +64,7 @@ struct efermi double ef_up = 0.0; ///< spin up Fermi energy double ef_dw = 0.0; ///< spin down Fermi energy bool two_efermi = false; ///< - double& get_ef(const int& is); + void set_efval(const int& is, const double& ef_in); double get_efval(const int& is) const; std::vector get_all_ef() const; }; diff --git a/source/module_elecstate/test/elecstate_fp_energy_test.cpp b/source/module_elecstate/test/elecstate_fp_energy_test.cpp index dc69b5a0bbf..36f775d9af9 100644 --- a/source/module_elecstate/test/elecstate_fp_energy_test.cpp +++ b/source/module_elecstate/test/elecstate_fp_energy_test.cpp @@ -14,7 +14,7 @@ * - fenergy::calculate_harris() * - fenergy::clear_all() * - fenergy::print_all() - * - efermi::get_ef() + * - efermi::set_efval() * - efermi::get_efval() */ class fenergy : public ::testing::Test @@ -61,19 +61,18 @@ TEST_F(fenergy, print_all) TEST_F(fenergy, eferm_get_ef) { eferm.two_efermi = false; - double& tmp_ef = eferm.get_ef(0); - tmp_ef = 0.7; + eferm.set_efval(0, 0.7); EXPECT_EQ(eferm.ef, 0.7); + eferm.set_efval(2, 0.77); + EXPECT_EQ(eferm.ef, 0.77); eferm.two_efermi = true; - double& tmp_efup = eferm.get_ef(0); - tmp_efup = 1.0; - EXPECT_EQ(eferm.ef_up, 1.0); - double& tmp_efdw = eferm.get_ef(1); - tmp_efdw = -1.0; + eferm.set_efval(0, 0.6); + EXPECT_EQ(eferm.ef_up, 0.6); + eferm.set_efval(1, -1.0); EXPECT_EQ(eferm.ef_dw, -1.0); testing::internal::CaptureStdout(); - EXPECT_EXIT(double& tmpp = eferm.get_ef(2);, ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(eferm.set_efval(3, 1.0);, ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("Please check NSPIN when TWO_EFERMI is true")); } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 8f3822cf529..68359163b5b 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -26,6 +26,7 @@ #include "module_io/nscf_band.h" #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" +#include "module_io/out_ldos.h" #include "module_io/to_wannier90_pw.h" #include "module_io/winput.h" #include "module_io/write_dos_pw.h" @@ -896,6 +897,11 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) } } + if (PARAM.inp.out_ldos) + { + ModuleIO::out_ldos(this->pelec, this->psi[0], this->Pgrid, ucell); + } + //! 4) Print out band structure information if (PARAM.inp.out_band[0]) { diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index c79e04e2eac..2166a15fb4e 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -290,7 +290,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) for (int is = 0; is < PARAM.inp.nspin; ++is) { - this->pelec->eferm.get_ef(is) = 0.; + this->pelec->eferm.set_efval(is, 0); this->theta_[is] = 0.; ModuleBase::GlobalFunc::ZEROS(this->pdLdphi_[is], this->pw_rho->nrxx); ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx); @@ -326,8 +326,7 @@ void ESolver_OF::update_potential(UnitCell& ucell) { this->pdEdphi_[is][ir] = vr_eff[ir]; } - this->pelec->eferm.get_ef(is) = this->cal_mu(this->pphi_[is], this->pdEdphi_[is], this->nelec_[is]); - + this->pelec->eferm.set_efval(is, this->cal_mu(this->pphi_[is], this->pdEdphi_[is], this->nelec_[is])); for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) { this->pdLdphi_[is][ir] @@ -540,8 +539,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso { this->pdEdphi_[0][ir] = vr_eff[ir]; } - this->pelec->eferm.get_ef(0) = this->cal_mu(this->pphi_[0], this->pdEdphi_[0], this->nelec_[0]); - + this->pelec->eferm.set_efval(0, this->cal_mu(this->pphi_[0], this->pdEdphi_[0], this->nelec_[0])); // === temporary === // assert(GlobalV::of_kinetic == "wt" || GlobalV::of_kinetic == "ml"); // ================= diff --git a/source/module_io/out_ldos.h b/source/module_io/out_ldos.h new file mode 100644 index 00000000000..4632338e6a0 --- /dev/null +++ b/source/module_io/out_ldos.h @@ -0,0 +1,27 @@ +#ifndef OUT_LDOS_H +#define OUT_LDOS_H + +#include "module_elecstate/elecstate_pw.h" +#include "module_io/cube_io.h" + +namespace ModuleIO +{ +void out_ldos(const elecstate::ElecState* pelec, + const psi::Psi>& psi, + const Parallel_Grid& pgrid, + const UnitCell& ucell) +{ + std::vector ldos(pelec->charge->nrxx); + + pelec->cal_ldos(psi, ldos); + + std::stringstream fn; + fn << PARAM.globalv.global_out_dir << "LDOS_" << PARAM.inp.stm_bias << "eV" + << ".cube"; + + ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, 11, 0); +} + +} // namespace ModuleIO + +#endif // OUT_LDOS_H From ce30e1b5457a1a286d4c7cddb15ffed857354dd2 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 2 Apr 2025 17:58:17 +0800 Subject: [PATCH 3/5] Refactor: after_all_runners --- source/module_elecstate/cal_ldos.cpp | 18 +++++- source/module_elecstate/elecstate_lcao.h | 3 + source/module_esolver/esolver.h | 2 +- source/module_esolver/esolver_fp.cpp | 8 +++ source/module_esolver/esolver_fp.h | 2 + source/module_esolver/esolver_ks.cpp | 69 +++++++++++++++++++++- source/module_esolver/esolver_ks.h | 2 + source/module_esolver/esolver_ks_lcao.cpp | 64 +-------------------- source/module_esolver/esolver_ks_pw.cpp | 70 ++--------------------- source/module_esolver/esolver_of.cpp | 6 +- source/module_io/out_ldos.h | 11 ++-- 11 files changed, 112 insertions(+), 143 deletions(-) diff --git a/source/module_elecstate/cal_ldos.cpp b/source/module_elecstate/cal_ldos.cpp index 8ef35a6f954..bf823f520ba 100644 --- a/source/module_elecstate/cal_ldos.cpp +++ b/source/module_elecstate/cal_ldos.cpp @@ -1,3 +1,4 @@ +#include "elecstate_lcao.h" #include "elecstate_pw.h" namespace elecstate @@ -9,7 +10,6 @@ void ElecStatePW::cal_ldos(const psi::Psi& psi, std::vecto const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0; const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0; - ldos.resize(this->basis->nrxx); std::vector wfcr(this->basis->nrxx); for (int ik = 0; ik < this->klist->get_nks(); ++ik) @@ -37,4 +37,20 @@ template class ElecStatePW, base_device::DEVICE_CPU>; template class ElecStatePW, base_device::DEVICE_GPU>; template class ElecStatePW, base_device::DEVICE_GPU>; #endif + +// lcao multi-k case +template <> +void ElecStateLCAO>::cal_ldos(const psi::Psi>& psi, + std::vector& ldos) const +{ +} + +// lcao Gamma_only case +template <> +void ElecStateLCAO::cal_ldos(const psi::Psi& psi, std::vector& ldos) const +{ +} + +template class ElecStateLCAO; // Gamma_only case +template class ElecStateLCAO>; // multi-k case } // namespace elecstate diff --git a/source/module_elecstate/elecstate_lcao.h b/source/module_elecstate/elecstate_lcao.h index 4beeb017f0f..9452b08426b 100644 --- a/source/module_elecstate/elecstate_lcao.h +++ b/source/module_elecstate/elecstate_lcao.h @@ -48,6 +48,9 @@ class ElecStateLCAO : public ElecState // const double* getRho(int spin) const override; virtual void cal_tau(const psi::Psi& psi) override; + //! calculate local density of states + virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const; + // update charge density for next scf step // void getNewRho() override; diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h index 800d40aeb9d..d8890523483 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -26,7 +26,7 @@ class ESolver virtual void runner(UnitCell& cell, const int istep) = 0; //! perform post processing calculations - virtual void after_all_runners(UnitCell& ucell){}; + virtual void after_all_runners(UnitCell& ucell) = 0; //! deal with exx and other calculation than scf/md/relax/cell-relax: //! such as nscf, get_wf and get_pchg diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index f99fb2264b7..78611955799 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -408,4 +408,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& } } +void ESolver_FP::after_all_runners(UnitCell& ucell) +{ + GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; + GlobalV::ofs_running << std::setprecision(16); + GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; + GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; +} + } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index a4507c7c553..3634c63be5a 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -49,6 +49,8 @@ class ESolver_FP: public ESolver //! Initialize of the first-principels energy solver virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override; + virtual void after_all_runners(UnitCell& ucell) override; + protected: //! Something to do before SCF iterations. virtual void before_scf(UnitCell& ucell, const int istep); diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index f06a0ce125d..4824af1b002 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -2,16 +2,18 @@ #include "module_base/timer.h" #include "module_cell/cal_atoms_info.h" +#include "module_elecstate/elecstate_print.h" #include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_hsolver/hsolver.h" #include "module_io/cube_io.h" #include "module_io/json_output/init_info.h" #include "module_io/json_output/output_info.h" +#include "module_io/nscf_band.h" +#include "module_io/out_ldos.h" #include "module_io/output_log.h" #include "module_io/print_info.h" #include "module_io/write_istate_info.h" #include "module_parameter/parameter.h" -#include "module_elecstate/elecstate_print.h" -#include "module_hsolver/hsolver.h" #include #include @@ -724,6 +726,69 @@ void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bo } } +template +void ESolver_KS::after_all_runners(UnitCell& ucell) +{ + ESolver_FP::after_all_runners(ucell); + + // 1) write information + if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0 || PARAM.inp.out_proj_band != 0) + { + GlobalV::ofs_running << "\n\n\n\n"; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>>>>>>>>>>>>>>" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " | Post-processing of data: " + " |" + << std::endl; + GlobalV::ofs_running << " | DOS (density of states) and bands will be " + "output here. |" + << std::endl; + GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken " + "charge analysis can be done. |" + << std::endl; + GlobalV::ofs_running << " | Also the .bxsf file containing fermi " + "surface information can be |" + << std::endl; + GlobalV::ofs_running << " | done here. " + " |" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" + "<<<<<<<<<<<<<<<<<<<<<<<<<" + << std::endl; + GlobalV::ofs_running << "\n\n\n\n"; + } + + // 2) write information + ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv); + + const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; + + // 3) print out band information + if (PARAM.inp.out_band[0]) + { + for (int is = 0; is < nspin0; is++) + { + std::stringstream ss2; + ss2 << PARAM.globalv.global_out_dir << "BANDS_" << is + 1 << ".dat"; + GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; + ModuleIO::nscf_band(is, + ss2.str(), + PARAM.inp.nbands, + 0.0, + PARAM.inp.out_band[1], + this->pelec->ekb, + this->kv); + } + } +} + //------------------------------------------------------------------------------ //! the 16th-20th functions of ESolver_KS //! mohan add 2024-05-12 diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index b0ef3eb279f..d1c2de4b17f 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -44,6 +44,8 @@ class ESolver_KS : public ESolver_FP virtual void runner(UnitCell& ucell, const int istep) override; + virtual void after_all_runners(UnitCell& ucell) override; + protected: //! Something to do before SCF iterations. virtual void before_scf(UnitCell& ucell, const int istep) override; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 941e4938be4..642b45b7742 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -12,7 +12,6 @@ #include "module_io/dos_nao.h" #include "module_io/io_dmk.h" #include "module_io/io_npz.h" -#include "module_io/nscf_band.h" #include "module_io/output_dmk.h" #include "module_io/output_log.h" #include "module_io/output_mat_sparse.h" @@ -393,71 +392,10 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); - GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; - GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; - GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; - - // 1) write information - if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0 || PARAM.inp.out_proj_band != 0) - { - GlobalV::ofs_running << "\n\n\n\n"; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" - ">>>>>>>>>>>>>>>>>>>>>>>>>" - << std::endl; - GlobalV::ofs_running << " | " - " |" - << std::endl; - GlobalV::ofs_running << " | Post-processing of data: " - " |" - << std::endl; - GlobalV::ofs_running << " | DOS (density of states) and bands will be " - "output here. |" - << std::endl; - GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken " - "charge analysis can be done. |" - << std::endl; - GlobalV::ofs_running << " | Also the .bxsf file containing fermi " - "surface information can be |" - << std::endl; - GlobalV::ofs_running << " | done here. " - " |" - << std::endl; - GlobalV::ofs_running << " | " - " |" - << std::endl; - GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" - "<<<<<<<<<<<<<<<<<<<<<<<<<" - << std::endl; - GlobalV::ofs_running << "\n\n\n\n"; - } - - // 2) write information - if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "md" || PARAM.inp.calculation == "relax") - { - ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv); - } + ESolver_KS::after_all_runners(ucell); const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; - // 3) print out band information - if (PARAM.inp.out_band[0]) - { - for (int is = 0; is < nspin0; is++) - { - std::stringstream ss2; - ss2 << PARAM.globalv.global_out_dir << "BANDS_" << is + 1 << ".dat"; - GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; - ModuleIO::nscf_band(is, - ss2.str(), - PARAM.inp.nbands, - 0.0, - PARAM.inp.out_band[1], - this->pelec->ekb, - this->kv); - } - } - // 4) write projected band structure by jiyy-2022-4-20 if (PARAM.inp.out_proj_band) { diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 68359163b5b..d5e7886b297 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -23,7 +23,6 @@ #include "module_hsolver/kernels/dngvd_op.h" #include "module_io/berryphase.h" #include "module_io/get_pchg_pw.h" -#include "module_io/nscf_band.h" #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" #include "module_io/out_ldos.h" @@ -826,53 +825,9 @@ void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& s template void ESolver_KS_PW::after_all_runners(UnitCell& ucell) { - //! 1) Output information to screen - GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; - GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; - GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; - - if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0) - { - GlobalV::ofs_running << "\n\n\n\n"; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" - ">>>>>>>>>>>>>>>>>>>>>>>>>" - << std::endl; - GlobalV::ofs_running << " | " - " |" - << std::endl; - GlobalV::ofs_running << " | Post-processing of data: " - " |" - << std::endl; - GlobalV::ofs_running << " | DOS (density of states) and bands will be " - "output here. |" - << std::endl; - GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken " - "charge analysis can be done. |" - << std::endl; - GlobalV::ofs_running << " | Also the .bxsf file containing fermi " - "surface information can be |" - << std::endl; - GlobalV::ofs_running << " | done here. " - " |" - << std::endl; - GlobalV::ofs_running << " | " - " |" - << std::endl; - GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" - "<<<<<<<<<<<<<<<<<<<<<<<<<" - << std::endl; - GlobalV::ofs_running << "\n\n\n\n"; - } - - int nspin0 = 1; - if (PARAM.inp.nspin == 2) - { - nspin0 = 2; - } - - //! 2) Print occupation numbers into istate.info - ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv); + ESolver_KS::after_all_runners(ucell); + + const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; //! 3) Compute density of states (DOS) if (PARAM.inp.out_dos) @@ -897,29 +852,12 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) } } + // out ldos if (PARAM.inp.out_ldos) { ModuleIO::out_ldos(this->pelec, this->psi[0], this->Pgrid, ucell); } - //! 4) Print out band structure information - if (PARAM.inp.out_band[0]) - { - for (int is = 0; is < nspin0; is++) - { - std::stringstream ss2; - ss2 << PARAM.globalv.global_out_dir << "BANDS_" << is + 1 << ".dat"; - GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; - ModuleIO::nscf_band(is, - ss2.str(), - PARAM.inp.nbands, - 0.0, - PARAM.inp.out_band[1], - this->pelec->ekb, - this->kv); - } - } - //! 5) Calculate the spillage value, used to generate numerical atomic orbitals if (PARAM.inp.basis_type == "pw" && winput::out_spillage) { diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 2166a15fb4e..c7e79d4729e 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -557,11 +557,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso */ void ESolver_OF::after_all_runners(UnitCell& ucell) { - - GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; - GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; - GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; + ESolver_FP::after_all_runners(ucell); } /** diff --git a/source/module_io/out_ldos.h b/source/module_io/out_ldos.h index 4632338e6a0..23e93c0f318 100644 --- a/source/module_io/out_ldos.h +++ b/source/module_io/out_ldos.h @@ -1,15 +1,16 @@ #ifndef OUT_LDOS_H #define OUT_LDOS_H -#include "module_elecstate/elecstate_pw.h" +#include "module_elecstate/elecstate.h" #include "module_io/cube_io.h" namespace ModuleIO { -void out_ldos(const elecstate::ElecState* pelec, - const psi::Psi>& psi, - const Parallel_Grid& pgrid, - const UnitCell& ucell) +template +inline void out_ldos(const elecstate::ElecState* pelec, + const psi::Psi& psi, + const Parallel_Grid& pgrid, + const UnitCell& ucell) { std::vector ldos(pelec->charge->nrxx); From eca763df0fc4a821c9158b241a153965cd788520 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 2 Apr 2025 18:52:02 +0800 Subject: [PATCH 4/5] Fix: complie bug --- source/module_elecstate/cal_ldos.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/source/module_elecstate/cal_ldos.cpp b/source/module_elecstate/cal_ldos.cpp index bf823f520ba..480a83d5de1 100644 --- a/source/module_elecstate/cal_ldos.cpp +++ b/source/module_elecstate/cal_ldos.cpp @@ -1,4 +1,6 @@ +#ifdef __LCAO #include "elecstate_lcao.h" +#endif #include "elecstate_pw.h" namespace elecstate @@ -38,6 +40,7 @@ template class ElecStatePW, base_device::DEVICE_GPU>; template class ElecStatePW, base_device::DEVICE_GPU>; #endif +#ifdef __LCAO // lcao multi-k case template <> void ElecStateLCAO>::cal_ldos(const psi::Psi>& psi, @@ -53,4 +56,5 @@ void ElecStateLCAO::cal_ldos(const psi::Psi& psi, std::vector; // Gamma_only case template class ElecStateLCAO>; // multi-k case +#endif } // namespace elecstate From 317c084c33e17afeb8ec45e8012f059f44b88c6b Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 2 Apr 2025 21:47:54 +0800 Subject: [PATCH 5/5] Refactor: cal_ldos in pw --- source/Makefile.Objects | 2 +- source/module_elecstate/CMakeLists.txt | 1 - source/module_elecstate/cal_ldos.cpp | 60 ------------------------ source/module_elecstate/elecstate.h | 9 ---- source/module_elecstate/elecstate_lcao.h | 3 -- source/module_elecstate/elecstate_pw.h | 7 +-- source/module_esolver/esolver_ks.cpp | 1 - source/module_esolver/esolver_ks_pw.cpp | 7 ++- source/module_io/CMakeLists.txt | 1 + source/module_io/cal_ldos.cpp | 55 ++++++++++++++++++++++ source/module_io/cal_ldos.h | 16 +++++++ source/module_io/out_ldos.h | 28 ----------- 12 files changed, 80 insertions(+), 110 deletions(-) delete mode 100644 source/module_elecstate/cal_ldos.cpp create mode 100644 source/module_io/cal_ldos.cpp create mode 100644 source/module_io/cal_ldos.h delete mode 100644 source/module_io/out_ldos.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index f99cf3613fd..8e258e8ebc5 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -235,7 +235,6 @@ OBJS_ELECSTAT=elecstate.o\ H_TDDFT_pw.o\ pot_xc.o\ cal_ux.o\ - cal_ldos.o\ read_orb.o\ cal_nelec_nband.o\ read_pseudo.o\ @@ -498,6 +497,7 @@ OBJS_IO=input_conv.o\ write_dos_pw.o\ nscf_band.o\ cal_dos.o\ + cal_ldos.o\ cif_io.o\ dos_nao.o\ numerical_descriptor.o\ diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 15e474789de..0bdd120879e 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -32,7 +32,6 @@ list(APPEND objects magnetism.cpp occupy.cpp cal_ux.cpp - cal_ldos.cpp read_orb.cpp cal_nelec_nband.cpp read_pseudo.cpp diff --git a/source/module_elecstate/cal_ldos.cpp b/source/module_elecstate/cal_ldos.cpp deleted file mode 100644 index 480a83d5de1..00000000000 --- a/source/module_elecstate/cal_ldos.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#ifdef __LCAO -#include "elecstate_lcao.h" -#endif -#include "elecstate_pw.h" - -namespace elecstate -{ -template -void ElecStatePW::cal_ldos(const psi::Psi& psi, std::vector& ldos) const -{ - // energy range for ldos (efermi as reference) - const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0; - const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0; - - std::vector wfcr(this->basis->nrxx); - - for (int ik = 0; ik < this->klist->get_nks(); ++ik) - { - psi.fix_k(ik); - double efermi = this->eferm.get_efval(this->klist->isk[ik]); - int nbands = psi.get_nbands(); - - for (int ib = 0; ib < nbands; ib++) - { - this->basis->recip2real(&psi(ib, 0), wfcr.data(), ik); - double eigenval = (ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV; - if (eigenval >= emin && eigenval <= emax) - { - for (int ir = 0; ir < basis->nrxx; ir++) - ldos[ir] += this->klist->wk[ik] * norm(wfcr[ir]); - } - } - } -} - -template class ElecStatePW, base_device::DEVICE_CPU>; -template class ElecStatePW, base_device::DEVICE_CPU>; -#if ((defined __CUDA) || (defined __ROCM)) -template class ElecStatePW, base_device::DEVICE_GPU>; -template class ElecStatePW, base_device::DEVICE_GPU>; -#endif - -#ifdef __LCAO -// lcao multi-k case -template <> -void ElecStateLCAO>::cal_ldos(const psi::Psi>& psi, - std::vector& ldos) const -{ -} - -// lcao Gamma_only case -template <> -void ElecStateLCAO::cal_ldos(const psi::Psi& psi, std::vector& ldos) const -{ -} - -template class ElecStateLCAO; // Gamma_only case -template class ElecStateLCAO>; // multi-k case -#endif -} // namespace elecstate diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index e83e714f381..5757ab696e5 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -61,15 +61,6 @@ class ElecState return; } - virtual void cal_ldos(const psi::Psi>& psi, std::vector& ldos) const - { - return; - } - virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const - { - return; - } - // update charge density for next scf step // in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge // density rho among these rho, diff --git a/source/module_elecstate/elecstate_lcao.h b/source/module_elecstate/elecstate_lcao.h index 9452b08426b..4beeb017f0f 100644 --- a/source/module_elecstate/elecstate_lcao.h +++ b/source/module_elecstate/elecstate_lcao.h @@ -48,9 +48,6 @@ class ElecStateLCAO : public ElecState // const double* getRho(int spin) const override; virtual void cal_tau(const psi::Psi& psi) override; - //! calculate local density of states - virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const; - // update charge density for next scf step // void getNewRho() override; diff --git a/source/module_elecstate/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h index da12f05a5c7..6007c8a4fe8 100644 --- a/source/module_elecstate/elecstate_pw.h +++ b/source/module_elecstate/elecstate_pw.h @@ -35,9 +35,6 @@ class ElecStatePW : public ElecState virtual void cal_tau(const psi::Psi& psi); - //! calculate local density of states - virtual void cal_ldos(const psi::Psi& psi, std::vector& ldos) const; - //! calculate becsum for uspp void cal_becsum(const psi::Psi& psi); @@ -49,12 +46,12 @@ class ElecStatePW : public ElecState T** rhog = nullptr; // [Device] [spin][nrxx] rhog Real** kin_r = nullptr; // [Device] [spin][nrxx] kin_r + ModulePW::PW_Basis_K* basis = nullptr; + protected: ModulePW::PW_Basis* rhopw_smooth = nullptr; - ModulePW::PW_Basis_K* basis = nullptr; - UnitCell* ucell = nullptr; const pseudopot_cell_vnl* ppcell = nullptr; diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 4824af1b002..4fa4ee0c022 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -9,7 +9,6 @@ #include "module_io/json_output/init_info.h" #include "module_io/json_output/output_info.h" #include "module_io/nscf_band.h" -#include "module_io/out_ldos.h" #include "module_io/output_log.h" #include "module_io/print_info.h" #include "module_io/write_istate_info.h" diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index d5e7886b297..23ffd297dd6 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -22,10 +22,10 @@ #include "module_hsolver/hsolver_pw.h" #include "module_hsolver/kernels/dngvd_op.h" #include "module_io/berryphase.h" +#include "module_io/cal_ldos.h" #include "module_io/get_pchg_pw.h" #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" -#include "module_io/out_ldos.h" #include "module_io/to_wannier90_pw.h" #include "module_io/winput.h" #include "module_io/write_dos_pw.h" @@ -855,7 +855,10 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) // out ldos if (PARAM.inp.out_ldos) { - ModuleIO::out_ldos(this->pelec, this->psi[0], this->Pgrid, ucell); + ModuleIO::cal_ldos(reinterpret_cast>*>(this->pelec), + this->psi[0], + this->Pgrid, + ucell); } //! 5) Calculate the spillage value, used to generate numerical atomic orbitals diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 677223d8444..7d8a396aba4 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -3,6 +3,7 @@ list(APPEND objects bessel_basis.cpp cal_test.cpp cal_dos.cpp + cal_ldos.cpp cif_io.cpp write_dos_pw.cpp nscf_band.cpp diff --git a/source/module_io/cal_ldos.cpp b/source/module_io/cal_ldos.cpp new file mode 100644 index 00000000000..7947ad1b470 --- /dev/null +++ b/source/module_io/cal_ldos.cpp @@ -0,0 +1,55 @@ +#include "cal_ldos.h" + +#include "cube_io.h" + +namespace ModuleIO +{ +void cal_ldos(const elecstate::ElecStatePW>* pelec, + const psi::Psi>& psi, + const Parallel_Grid& pgrid, + const UnitCell& ucell) +{ + // energy range for ldos (efermi as reference) + const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0; + const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0; + + std::vector ldos(pelec->charge->nrxx); + std::vector> wfcr(pelec->basis->nrxx); + + for (int ik = 0; ik < pelec->klist->get_nks(); ++ik) + { + psi.fix_k(ik); + double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]); + int nbands = psi.get_nbands(); + + for (int ib = 0; ib < nbands; ib++) + { + pelec->basis->recip2real(&psi(ib, 0), wfcr.data(), ik); + double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV; + if (eigenval >= emin && eigenval <= emax) + { + for (int ir = 0; ir < pelec->basis->nrxx; ir++) + ldos[ir] += pelec->klist->wk[ik] * norm(wfcr[ir]); + } + } + } + + std::stringstream fn; + fn << PARAM.globalv.global_out_dir << "LDOS_" << PARAM.inp.stm_bias << "eV" + << ".cube"; + + ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, 11, 0); +} + +#ifdef __LCAO +// lcao multi-k case +// void cal_ldos(elecstate::ElecState* pelec, const psi::Psi>& psi, std::vector& ldos) +// { +// } + +// // lcao Gamma_only case +// void cal_ldos(elecstate::ElecState* pelec, const psi::Psi& psi, std::vector& ldos) +// { +// } +#endif +} // namespace elecstate diff --git a/source/module_io/cal_ldos.h b/source/module_io/cal_ldos.h new file mode 100644 index 00000000000..e761b6f2002 --- /dev/null +++ b/source/module_io/cal_ldos.h @@ -0,0 +1,16 @@ +#ifndef CAL_LDOS_H +#define CAL_LDOS_H + +#include "module_elecstate/elecstate_pw.h" + +namespace ModuleIO +{ + +void cal_ldos(const elecstate::ElecStatePW>* pelec, + const psi::Psi>& psi, + const Parallel_Grid& pgrid, + const UnitCell& ucell); + +} // namespace ModuleIO + +#endif // CAL_LDOS_H \ No newline at end of file diff --git a/source/module_io/out_ldos.h b/source/module_io/out_ldos.h deleted file mode 100644 index 23e93c0f318..00000000000 --- a/source/module_io/out_ldos.h +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef OUT_LDOS_H -#define OUT_LDOS_H - -#include "module_elecstate/elecstate.h" -#include "module_io/cube_io.h" - -namespace ModuleIO -{ -template -inline void out_ldos(const elecstate::ElecState* pelec, - const psi::Psi& psi, - const Parallel_Grid& pgrid, - const UnitCell& ucell) -{ - std::vector ldos(pelec->charge->nrxx); - - pelec->cal_ldos(psi, ldos); - - std::stringstream fn; - fn << PARAM.globalv.global_out_dir << "LDOS_" << PARAM.inp.stm_bias << "eV" - << ".cube"; - - ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, 11, 0); -} - -} // namespace ModuleIO - -#endif // OUT_LDOS_H