Feature: ldos in pw basis set (#6100)

* add input parameters for ldos

* Feature: ldos in pw basis set

* Refactor: after_all_runners

* Fix: complie bug

* Refactor: cal_ldos in pw

Author: YuLiu98

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)
@@ -175,6 +176,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)
@@ -1707,6 +1709,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)
@@ -1963,9 +1971,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

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\

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;

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
     {

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<double> get_all_ef() const;
 };

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"));
 }

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

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

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);

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 <ctime>
 #include <iostream>
@@ -724,6 +725,69 @@ void ESolver_KS<T, Device>::after_scf(UnitCell& ucell, const int istep, const bo
     }
 }
 
+template <typename T, typename Device>
+void ESolver_KS<T, Device>::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