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/Makefile.Objects b/source/Makefile.Objects index 9b690c5fdf8..8e258e8ebc5 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -497,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/elecstate_pw.h b/source/module_elecstate/elecstate_pw.h index 679b9b712cb..6007c8a4fe8 100644 --- a/source/module_elecstate/elecstate_pw.h +++ b/source/module_elecstate/elecstate_pw.h @@ -46,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_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.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..4fa4ee0c022 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -2,16 +2,17 @@ #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/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 +725,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 8f3822cf529..23ffd297dd6 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -22,8 +22,8 @@ #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/nscf_band.h" #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" #include "module_io/to_wannier90_pw.h" @@ -825,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) @@ -896,22 +852,13 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) } } - //! 4) Print out band structure information - if (PARAM.inp.out_band[0]) + // out ldos + if (PARAM.inp.out_ldos) { - 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); - } + 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_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index c79e04e2eac..c7e79d4729e 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"); // ================= @@ -559,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/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/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