Merge branch 'develop' into fft12

Author: A-006

docs/advanced/acceleration/cuda.md

@@ -46,6 +46,7 @@ We provides [examples](https://github.com/deepmodeling/abacus-develop/tree/devel
 PW basis:
 - Only k point parallelization is supported, so the input keyword `kpar` will be set to match the number of MPI tasks automatically.
 - By default, CUDA architectures 60, 70, 75, 80, 86, and 89 are compiled (if supported). It can be overriden using the CMake variable [`CMAKE_CUDA_ARCHITECTURES`](https://cmake.org/cmake/help/latest/variable/CMAKE_CUDA_ARCHITECTURES.html) or the environmental variable [`CUDAARCHS`](https://cmake.org/cmake/help/latest/envvar/CUDAARCHS.html).
+
 LCAO basis:
 - Unless there is a specific reason, avoid using multiple GPUs, as it can be slower than using a single GPU. This is because the generalized eigenvalue solution of the LCAO basis set will incur additional communication overhead when calculated on multiple cards. When the memory limit of a GPU card makes it insufficient to complete the task, it is recommended to use multiple cards for calculation.
 - When using elpa on GPUs, some ELPA internal logs will be output.

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)
@@ -154,6 +155,7 @@
     - [out\_mat\_dh](#out_mat_dh)
     - [out\_mat\_xc](#out_mat_xc)
     - [out\_mat\_xc2](#out_mat_xc2)
+    - [out\_mat\_l](#out_mat_l)
     - [out\_eband\_terms](#out_eband_terms)
     - [out\_hr\_npz/out\_dm\_npz](#out_hr_npzout_dm_npz)
     - [dm\_to\_rho](#dm_to_rho)
@@ -174,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)
@@ -1706,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)
@@ -1807,6 +1816,13 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 - **Description**: Whether to print the exchange-correlation matrices in **numerical orbital representation** (unit: Ry): $\braket{\phi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\phi_j}(\mathbf{R})$ in CSR format (the same format as [out_mat_hs2](../elec_properties/hs_matrix.md#out_mat_hs2)) in the directory `OUT.${suffix}`. (Note that currently DeePKS term is not included. ) The files are named `Vxc_R_spin$s`.
 - **Default**: False
 
+### out_mat_l
+
+- **Type**: Boolean [Integer\](optional)
+- **Availability**: Numerical atomic orbital (NAO) basis
+- **Description**: Whether to print the expectation value of the angular momentum operator $\hat{L}_x$, $\hat{L}_y$, and $\hat{L}_z$ in the basis of the localized atomic orbitals. The files are named `OUT.${suffix}/${suffix}_Lx.dat`, `OUT.${suffix}/${suffix}_Ly.dat`, and `OUT.${suffix}/${suffix}_Lz.dat`. The second integer controls the precision of the output.
+- **Default**: False 8
+
 ### out_eband_terms
 
 - **Type**: Boolean
@@ -1955,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\
@@ -554,6 +555,7 @@ OBJS_IO=input_conv.o\
     read_input_item_output.o\
     read_set_globalv.o\
     orb_io.o\
+    cal_pLpR.o\
 
 OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_orb_info.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);