Merge branch 'develop' into pulay_fs

.github/workflows/coverage.yml

@@ -15,10 +15,10 @@ jobs:
         uses: actions/checkout@v4
       - name: Install Requirements for Coverage Testing
         run: |
-          apt update && apt install -y lcov
+          apt update && apt install -y lcov gpg
       - name: Building
         run: |
-          cmake -B build -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DBUILD_TESTING=ON -DENABLE_COVERAGE=ON
+          cmake -B build -DENABLE_COVERAGE=ON -DBUILD_TESTING=ON -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_PAW=ON -DENABLE_GOOGLEBENCH=ON -DENABLE_RAPIDJSON=ON
           cmake --build build -j`nproc`
           cmake --install build
       - name: Testing

.github/workflows/test.yml

@@ -39,7 +39,7 @@ jobs:
             --from-ref ${{ github.event.pull_request.base.sha }}
               --to-ref ${{ github.event.pull_request.head.sha }}
         continue-on-error: true
-      - uses: pre-commit-ci/lite-action@v1.0.3
+      - uses: pre-commit-ci/lite-action@v1.1.0
 
       - name: Build
         run: |

docs/advanced/install.md

@@ -93,9 +93,9 @@ cmake -B build -DUSE_CUDA=1 -DCMAKE_CUDA_COMPILER=${path to cuda toolkit}/bin/nv
 
 ## Build math library from source
 
-> Note: This flag is **enabled by default**. It will get better performance than the standard implementation on `gcc` and `clang`. But it **will be disabled** when using `Intel Compiler` since the math functions will get wrong results and the performance is also unexpectly poor.
+> Note: We recommend using the latest available compiler sets, since they offer faster implementations of math functions.
 
-To build math functions from source code, instead of using c++ standard implementation, define `USE_ABACUS_LIBM` flag.
+This flag is disabled by default. To build math functions from source code, define `USE_ABACUS_LIBM` flag. It is expected to get a better performance on legacy versions of `gcc` and `clang`.
 
 Currently supported math functions:
  `sin`, `cos`, `sincos`, `exp`, `cexp`

source/driver.cpp

@@ -41,7 +41,7 @@ void Driver::init()
 
     // (3) output information
     time_t time_finish = std::time(nullptr);
-    Print_Info::print_time(time_start, time_finish);
+    ModuleIO::print_time(time_start, time_finish);
 
     // (4) close all of the running logs
     ModuleBase::Global_File::close_all_log(GlobalV::MY_RANK, PARAM.inp.out_alllog,PARAM.inp.calculation);

source/module_elecstate/elecstate_lcao_tddft.cpp

@@ -22,13 +22,8 @@ void ElecStateLCAO_TDDFT::psiToRho_td(const psi::Psi<std::complex<double>>& psi)
 
     // this part for calculating DMK in 2d-block format, not used for charge now
     //    psi::Psi<std::complex<double>> dm_k_2d();
-
-    if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "scalapack_gvx"
-        || PARAM.inp.ks_solver == "lapack") // Peize Lin test 2019-05-15
-    {
-        elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM));
-        this->DM->cal_DMR();
-    }
+    elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM));
+    this->DM->cal_DMR();
 
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -227,7 +227,7 @@ int H_TDDFT_pw::check_ncut(int t_type)
         break;
 
     case 3:
-        ncut = 1;
+        ncut = 2;
         break;
 
         // case 4:

source/module_esolver/esolver.h

@@ -58,11 +58,8 @@ class ESolver
         return 0;
     }
 
-    // get conv_elec used in current scf
-    virtual bool get_conv_elec()
-    {
-        return true;
-    }
+    bool conv_esolver = true; // whether esolver is converged
+
     std::string classname;
 };
 

source/module_esolver/esolver_fp.cpp

@@ -119,23 +119,14 @@ void ESolver_FP::before_all_runners(const Input_para& inp, UnitCell& cell)
     return;
 }
 
-//------------------------------------------------------------------------------
-//! the 12th function of ESolver_KS: get_conv_elec
-//! tqzhao add 2024-05-15
-//------------------------------------------------------------------------------
-bool ESolver_FP::get_conv_elec()
-{
-    return this->conv_elec;
-}
-
 //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
 void ESolver_FP::after_scf(const int istep)
 {
     // 0) output convergence information
-    ModuleIO::output_convergence_after_scf(this->conv_elec, this->pelec->f_en.etot);
+    ModuleIO::output_convergence_after_scf(this->conv_esolver, this->pelec->f_en.etot);
 
     // 1) write fermi energy
-    ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef);
+    ModuleIO::output_efermi(this->conv_esolver, this->pelec->eferm.ef);
 
     // 2) update delta rho for charge extrapolation
     CE.update_delta_rho(GlobalC::ucell, &(this->chr), &(this->sf));

source/module_esolver/esolver_fp.h

@@ -44,9 +44,6 @@ namespace ModuleESolver
         //! Initialize of the first-principels energy solver
         virtual void before_all_runners(const Input_para& inp, UnitCell& cell) override;
 
-        // get conv_elec used in current scf
-        virtual bool get_conv_elec() override;
-
         virtual void init_after_vc(const Input_para& inp, UnitCell& cell);    // liuyu add 2023-03-09
 
         //! Electronic states
@@ -64,8 +61,6 @@ namespace ModuleESolver
         //! K points in Brillouin zone
         K_Vectors kv;
 
-        bool conv_elec; // If electron density is converged in scf.
-
       protected:
         //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
         virtual void after_scf(const int istep);

source/module_esolver/esolver_ks.cpp

@@ -7,11 +7,12 @@
 #include <chrono>
 #endif
 #include "module_base/timer.h"
+#include "module_cell/cal_atoms_info.h"
 #include "module_io/json_output/init_info.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_cell/cal_atoms_info.h"
 
 #include <iostream>
 //--------------Temporary----------------
@@ -221,7 +222,7 @@ void ESolver_KS<T, Device>::before_all_runners(const Input_para& inp, UnitCell&
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
     //! 7) print information
-    Print_Info::setup_parameters(ucell, this->kv);
+    ModuleIO::setup_parameters(ucell, this->kv);
 
     //! 8) new plane wave basis, fft grids, etc.
 #ifdef __MPI
@@ -419,7 +420,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
     bool firstscf = true;
-    this->conv_elec = false;
+    this->conv_esolver = false;
     this->niter = this->maxniter;
 
     // 4) SCF iterations
@@ -429,7 +430,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
     for (int iter = 1; iter <= this->maxniter; ++iter)
     {
         // 5) write head
-        this->write_head(GlobalV::ofs_running, istep, iter);
+        ModuleIO::write_head(GlobalV::ofs_running, istep, iter, this->basisname);
 
 #ifdef __MPI
         auto iterstart = MPI_Wtime();
@@ -496,9 +497,9 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
                                                           PARAM.inp.esolver_type,
                                                           diag_ethr,
                                                           PARAM.inp.nelec);
-                
+
                 // The error of HSolver is larger than drho,
-                // so a more precise HSolver should be excuconv_elected.
+                // so a more precise HSolver should be executed.
                 if (hsolver_error > drho)
                 {
                     diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running,
@@ -540,11 +541,11 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
             }
 #endif
 
-            this->conv_elec = (drho < this->scf_thr && not_restart_step && is_U_converged);
+            this->conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged);
 
             // If drho < hsolver_error in the first iter or drho < scf_thr, we
             // do not change rho.
-            if (drho < hsolver_error || this->conv_elec)
+            if (drho < hsolver_error || this->conv_esolver)
             {
                 if (drho < hsolver_error)
                 {
@@ -577,13 +578,13 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
         }
 #ifdef __MPI
         MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, PARAPW_WORLD);
-        MPI_Bcast(&this->conv_elec, 1, MPI_DOUBLE, 0, PARAPW_WORLD);
+        MPI_Bcast(&this->conv_esolver, 1, MPI_DOUBLE, 0, PARAPW_WORLD);
         MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, PARAPW_WORLD);
 #endif
 
         // 9) update potential
         // Hamilt should be used after it is constructed.
-        // this->phamilt->update(conv_elec);
+        // this->phamilt->update(conv_esolver);
         this->update_pot(istep, iter);
 
         // 10) finish scf iterations
@@ -617,7 +618,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
 #endif //__RAPIDJSON
 
         // 13) check convergence
-        if (this->conv_elec)
+        if (this->conv_esolver)
         {
             this->niter = iter;
             break;
@@ -632,7 +633,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
     std::cout << " >> Leave SCF iteration.\n * * * * * *" << std::endl;
 #ifdef __RAPIDJSON
     // 14) add Json of efermi energy converge
-    Json::add_output_efermi_converge(this->pelec->eferm.ef * ModuleBase::Ry_to_eV, this->conv_elec);
+    Json::add_output_efermi_converge(this->pelec->eferm.ef * ModuleBase::Ry_to_eV, this->conv_esolver);
 #endif //__RAPIDJSON
 
     // 15) after scf
@@ -666,9 +667,9 @@ void ESolver_KS<T, Device>::iter_finish(int& iter)
     this->pelec->f_en.etot_old = this->pelec->f_en.etot;
 
     // add a energy threshold for SCF convergence
-    if (this->conv_elec == 0) // only check when density is not converged
+    if (this->conv_esolver == 0) // only check when density is not converged
     {
-        this->conv_elec
+        this->conv_esolver
             = (iter != 1 && std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr);
     }
 }
@@ -687,33 +688,6 @@ void ESolver_KS<T, Device>::after_scf(const int istep)
     }
 }
 
-//------------------------------------------------------------------------------
-//! the 8th function of ESolver_KS: print_head
-//! mohan add 2024-05-12
-//------------------------------------------------------------------------------
-template <typename T, typename Device>
-void ESolver_KS<T, Device>::print_head()
-{
-    std::cout << " " << std::setw(7) << "ITER";
-
-    if (PARAM.inp.nspin == 2)
-    {
-        std::cout << std::setw(10) << "TMAG";
-        std::cout << std::setw(10) << "AMAG";
-    }
-
-    std::cout << std::setw(15) << "ETOT(eV)";
-    std::cout << std::setw(15) << "EDIFF(eV)";
-    std::cout << std::setw(11) << "DRHO";
-
-    if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
-    {
-        std::cout << std::setw(11) << "DKIN";
-    }
-
-    std::cout << std::setw(11) << "TIME(s)" << std::endl;
-}
-
 //------------------------------------------------------------------------------
 //! the 8th function of ESolver_KS: print_iter
 //! mohan add 2024-05-12
@@ -725,18 +699,7 @@ void ESolver_KS<T, Device>::print_iter(const int iter,
                                        const double duration,
                                        const double ethr)
 {
-    this->pelec->print_etot(this->conv_elec, iter, drho, dkin, duration, PARAM.inp.printe, ethr);
-}
-
-//------------------------------------------------------------------------------
-//! the 9th function of ESolver_KS: write_head
-//! mohan add 2024-05-12
-//------------------------------------------------------------------------------
-template <typename T, typename Device>
-void ESolver_KS<T, Device>::write_head(std::ofstream& ofs_running, const int istep, const int iter)
-{
-    ofs_running << "\n " << this->basisname << " ALGORITHM --------------- ION=" << std::setw(4) << istep + 1
-                << "  ELEC=" << std::setw(4) << iter << "--------------------------------\n";
+    this->pelec->print_etot(this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, ethr);
 }
 
 //------------------------------------------------------------------------------

source/module_esolver/esolver_ks.h

@@ -73,12 +73,7 @@ class ESolver_KS : public ESolver_FP
 		virtual void update_pot(const int istep, const int iter) {};
 
     protected:
-
-		// Print the headline on the screen:
-		// ITER   ETOT(eV)       EDIFF(eV)      DRHO    TIME(s) 
-		void print_head();
-
-		// Print inforamtion in each iter
+        // Print inforamtion in each iter
 		// G1    -3.435545e+03  0.000000e+00   3.607e-01  2.862e-01
 		// for metaGGA
 		// ITER   ETOT(eV)       EDIFF(eV)      DRHO       DKIN       TIME(s) 
@@ -90,13 +85,6 @@ class ESolver_KS : public ESolver_FP
 				const double duration, 
 				const double ethr);
 
-		// Write the headline in the running_log file
-		// "PW/LCAO" ALGORITHM --------------- ION=   1  ELEC=   1--------------------------------
-		void write_head(
-				std::ofstream& ofs_running, 
-				const int istep, 
-				const int iter);
-
         //! Hamiltonian
 		hamilt::Hamilt<T, Device>* p_hamilt = nullptr;
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -136,7 +136,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
         this->kv.set(ucell.symm, PARAM.inp.kpoint_file, PARAM.inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
-        Print_Info::setup_parameters(ucell, this->kv);
+        ModuleIO::setup_parameters(ucell, this->kv);
     }
     else
     {
@@ -817,7 +817,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
     ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot");
 
     // 1) print Hamiltonian and Overlap matrix
-    if (this->conv_elec || iter == PARAM.inp.scf_nmax)
+    if (this->conv_esolver || iter == PARAM.inp.scf_nmax)
     {
         if (!PARAM.globalv.gamma_only_local && (PARAM.inp.out_mat_hs[0] || PARAM.inp.deepks_v_delta))
         {
@@ -875,7 +875,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
     }
 
     // 2) print wavefunctions
-    if (elecstate::ElecStateLCAO<TK>::out_wfc_lcao && (this->conv_elec || iter == PARAM.inp.scf_nmax)
+    if (elecstate::ElecStateLCAO<TK>::out_wfc_lcao && (this->conv_esolver || iter == PARAM.inp.scf_nmax)
         && (istep % PARAM.inp.out_interval == 0))
     {
         ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO<TK>::out_wfc_lcao,
@@ -887,7 +887,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
                                 istep);
     }
 
-    if (!this->conv_elec)
+    if (!this->conv_esolver)
     {
         if (PARAM.inp.nspin == 4)
         {
@@ -977,7 +977,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
         }
     }
 
-    if (GlobalC::exx_info.info_global.cal_exx && this->conv_elec)
+    if (GlobalC::exx_info.info_global.cal_exx && this->conv_esolver)
     {
         // Kerker mixing does not work for the density matrix.
         // In the separate loop case, it can still work in the subsequent inner loops where Hexx(DM) is fixed.
@@ -989,7 +989,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
         }
         if (GlobalC::exx_info.info_ri.real_number)
         {
-            this->conv_elec = this->exd->exx_after_converge(
+            this->conv_esolver = this->exd->exx_after_converge(
                 *this->p_hamilt,
                 *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                 this->kv,
@@ -1000,7 +1000,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
         }
         else
         {
-            this->conv_elec = this->exc->exx_after_converge(
+            this->conv_esolver = this->exc->exx_after_converge(
                 *this->p_hamilt,
                 *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                 this->kv,
@@ -1080,7 +1080,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
     }
 
     // 6) use the converged occupation matrix for next MD/Relax SCF calculation
-    if (PARAM.inp.dft_plus_u && this->conv_elec)
+    if (PARAM.inp.dft_plus_u && this->conv_esolver)
     {
         GlobalC::dftu.initialed_locale = true;
     }