Test: add test for pw LDOS

source/module_esolver/esolver_ks_lcao.cpp

@@ -8,6 +8,7 @@
 #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_io/berryphase.h"
+#include "module_io/cal_ldos.h"
 #include "module_io/cube_io.h"
 #include "module_io/dos_nao.h"
 #include "module_io/io_dmk.h"
@@ -419,6 +420,15 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                               this->p_hamilt);
     }
 
+    // out ldos
+    if (PARAM.inp.out_ldos[0])
+    {
+        ModuleIO::Cal_ldos<TK>::cal_ldos_lcao(reinterpret_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec),
+                                              this->psi[0],
+                                              this->Pgrid,
+                                              ucell);
+    }
+
     // 6) print out exchange-correlation potential
     if (PARAM.inp.out_mat_xc)
     {

source/module_esolver/esolver_ks_pw.cpp

@@ -853,12 +853,13 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
     }
 
     // out ldos
-    if (PARAM.inp.out_ldos)
+    if (PARAM.inp.out_ldos[0])
     {
-        ModuleIO::cal_ldos(reinterpret_cast<elecstate::ElecStatePW<std::complex<double>>*>(this->pelec),
-                           this->psi[0],
-                           this->Pgrid,
-                           ucell);
+        ModuleIO::Cal_ldos<std::complex<double>>::cal_ldos_pw(
+            reinterpret_cast<elecstate::ElecStatePW<std::complex<double>>*>(this->pelec),
+            this->psi[0],
+            this->Pgrid,
+            ucell);
     }
 
     //! 5) Calculate the spillage value, used to generate numerical atomic orbitals

source/module_io/cal_ldos.cpp

@@ -1,13 +1,18 @@
 #include "cal_ldos.h"
 
 #include "cube_io.h"
+#include "module_base/blas_connector.h"
+#include "module_base/scalapack_connector.h"
+
+#include <type_traits>
 
 namespace ModuleIO
 {
-void cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
-              const psi::Psi<std::complex<double>>& psi,
-              const Parallel_Grid& pgrid,
-              const UnitCell& ucell)
+template <typename T>
+void Cal_ldos<T>::cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>* pelec,
+                              const psi::Psi<std::complex<double>>& 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;
@@ -19,13 +24,13 @@ void cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
     for (int ik = 0; ik < pelec->klist->get_nks(); ++ik)
     {
         psi.fix_k(ik);
-        double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
+        const 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;
+            const 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++)
@@ -38,18 +43,90 @@ void cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
     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);
+    const int precision = PARAM.inp.out_ldos[1];
+    ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0);
 }
 
 #ifdef __LCAO
-// lcao multi-k case
-// void cal_ldos(elecstate::ElecState* pelec, const psi::Psi<std::complex<double>>& psi, std::vector<double>& ldos)
-// {
-// }
-
-// // lcao Gamma_only case
-// void cal_ldos(elecstate::ElecState* pelec, const psi::Psi<double>& psi, std::vector<double>& ldos)
-// {
-// }
+template <typename T>
+void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
+                                const psi::Psi<T>& 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;
+
+    // calulate dm-like
+    const int nbands_local = psi.get_nbands();
+    const int nbasis_local = psi.get_nbasis();
+
+    // psi.T * wk * psi.conj()
+    // result[ik](iw1,iw2) = \sum_{ib} psi[ik](ib,iw1).T * wk(k) * psi[ik](ib,iw2).conj()
+    for (int ik = 0; ik < psi.get_nk(); ++ik)
+    {
+        psi.fix_k(ik);
+        const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
+
+        // T* dmk_pointer = DM.get_DMK_pointer(ik);
+
+        psi::Psi<T> wk_psi(1, psi.get_nbands(), psi.get_nbasis(), psi.get_nbasis(), true);
+        const T* ppsi = psi.get_pointer();
+        T* pwk_psi = wk_psi.get_pointer();
+
+        // #ifdef _OPENMP
+        // #pragma omp parallel for schedule(static, 1024)
+        // #endif
+        //         for (int i = 0; i < wk_psi.size(); ++i)
+        //         {
+        //             pwk_psi[i] = my_conj(ppsi[i]);
+        //         }
+
+        //         int ib_global = 0;
+        //         for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
+        //         {
+        //             while (ib_local != ParaV->global2local_col(ib_global))
+        //             {
+        //                 ++ib_global;
+        //                 if (ib_global >= wg.nc)
+        //                 {
+        //                     ModuleBase::WARNING_QUIT("cal_ldos", "please check global2local_col!");
+        //                 }
+        //             }
+
+        //             const double eigenval = (pelec->ekb(ik, ib_global) - 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]);
+        //             }
+
+        //             double* wg_wfc_pointer = &(wk_psi(0, ib_local, 0));
+        //             BlasConnector::scal(nbasis_local, pelec->klist->wk[ik], wg_wfc_pointer, 1);
+        //         }
+
+        //         // C++: dm(iw1,iw2) = psi(ib,iw1).T * wk_psi(ib,iw2)
+        // #ifdef __MPI
+        //         psiMulPsiMpi(wk_psi, psi, dmk_pointer, ParaV->desc_wfc, ParaV->desc);
+        // #else
+        //         psiMulPsi(wk_psi, psi, dmk_pointer);
+        // #endif
+    }
+}
+
+double my_conj(double x)
+{
+    return x;
+}
+
+std::complex<double> my_conj(const std::complex<double>& z)
+{
+    return {z.real(), -z.imag()};
+}
+
 #endif
+
+template class Cal_ldos<double>;               // Gamma_only case
+template class Cal_ldos<std::complex<double>>; // multi-k case
 } // namespace elecstate

source/module_io/cal_ldos.h

@@ -1,16 +1,29 @@
 #ifndef CAL_LDOS_H
 #define CAL_LDOS_H
 
+#include "module_elecstate/elecstate_lcao.h"
 #include "module_elecstate/elecstate_pw.h"
 
 namespace ModuleIO
 {
+template <typename T>
+class Cal_ldos
+{
+  public:
+    Cal_ldos(){};
+    ~Cal_ldos(){};
+
+    static void cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>* pelec,
+                            const psi::Psi<std::complex<double>>& psi,
+                            const Parallel_Grid& pgrid,
+                            const UnitCell& ucell);
 
-void cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
-              const psi::Psi<std::complex<double>>& psi,
-              const Parallel_Grid& pgrid,
-              const UnitCell& ucell);
+    static void cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
+                              const psi::Psi<T>& psi,
+                              const Parallel_Grid& pgrid,
+                              const UnitCell& ucell);
 
+}; // namespace Cal_ldos
 } // namespace ModuleIO
 
 #endif // CAL_LDOS_H

source/module_io/read_input_item_output.cpp

@@ -145,8 +145,17 @@ void ReadInput::item_output()
     }
     {
         Input_Item item("out_ldos");
-        item.annotation = "output local density of states";
-        read_sync_bool(input.out_ldos);
+        item.annotation = "output local density of states, second parameter controls the precision";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            const size_t count = item.get_size();
+            if (count != 1 && count != 2)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "out_ldos should have 1 or 2 values");
+            }
+            para.input.out_ldos[0] = assume_as_boolean(item.str_values[0]);
+            para.input.out_ldos[1] = (count == 2) ? std::stoi(item.str_values[1]) : 3;
+        };
+        sync_intvec(input.out_ldos, 2, 0);
         this->add_item(item);
     }
     {

source/module_io/read_input_item_postprocess.cpp

@@ -58,7 +58,7 @@ void ReadInput::item_postprocess()
         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)
+            if (para.input.out_ldos[0] && para.input.stm_bias == 0.0)
             {
                 ModuleBase::WARNING_QUIT("ReadInput", "a nonzero stm_bias is required for ldos calculation");
             }

source/module_io/test/read_input_ptest.cpp

@@ -201,7 +201,8 @@ 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_ldos[0], 1);
+    EXPECT_EQ(param.inp.out_ldos[1], 3);
     EXPECT_EQ(param.inp.out_band[0], 0);
     EXPECT_EQ(param.inp.out_band[1], 8);
     EXPECT_EQ(param.inp.out_proj_band, 0);

source/module_io/test/support/INPUT

@@ -65,7 +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_ldos                       1 #output local density of states, second parameter controls the precision
 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

source/module_io/test_serial/read_input_item_test.cpp

@@ -393,7 +393,7 @@ TEST_F(InputTest, Item_test)
     }
     { // stm_bias
         auto it = find_label("stm_bias", readinput.input_lists);
-        param.input.out_ldos = true;
+        param.input.out_ldos[0] = 1;
         param.input.stm_bias = 0.0;
 
         testing::internal::CaptureStdout();

source/module_parameter/input_parameter.h

@@ -363,7 +363,7 @@ struct Input_para
     int printe = 0;                       ///< Print out energy for each band for every printe step, default is scf_nmax
     std::vector<int> 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
+    std::vector<int> out_ldos = {0, 3};   ///< 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.

tests/integrate/101_PW_15_pseudopots/INPUT

@@ -23,3 +23,6 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+
+out_ldos              1
+stm_bias              2