Merge branch 'develop' of https://github.com/mohanchen/abacus-mc into develop

Author: mohanchen

docs/advanced/input_files/input-main.md

@@ -599,7 +599,7 @@ These variables are used to control general system parameters.
   - nao: from numerical atomic orbitals. If they are not enough, other wave functions are initialized with random numbers.
   - nao+random: add small random numbers on numerical atomic orbitals
 
-  > Only the `file` option is useful for the lcao basis set, which is mostly used when [calculation](#calculation) is set to `set_wf` and `get_pchg`. See more details in [out_wfc_lcao](#out_wfc_lcao).
+  > Only the `file` option is useful for the lcao basis set, which is mostly used when [calculation](#calculation) is set to `get_wf` and `get_pchg`. See more details in [out_wfc_lcao](#out_wfc_lcao).
 - **Default**: atomic
 
 ### init_chg
@@ -1928,7 +1928,7 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 ### if_separate_k
 
 - **Type**: Boolean
-- **Availability**: Only for LCAO, used only when `calculation = get_pchg` and `gamma_only` is turned off.
+- **Availability**: For both PW and LCAO. When `basis_type = pw`, used if `out_pchg` is set. When `basis_type = lcao`, used only when `calculation = get_pchg` and `gamma_only` is turned off.
 - **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. **Warning**: Enabling symmetry may produce incorrect results due to incorrect k-point weights. Therefore, when calculating partial charge densities, it is strongly recommended to set `symmetry = -1`.
 - **Default**: false
 

source/module_elecstate/module_charge/symmetry_rho.cpp

@@ -97,7 +97,7 @@ void Symmetry_rho::psymm(double* rho_part,
         rhotot.resize(rho_basis->nxyz);
         ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz);
     }
-    Pgrid.reduce(rhotot.data(), rho_part);
+    Pgrid.reduce(rhotot.data(), rho_part, false);
 
     // (2)
     if (GlobalV::MY_RANK == 0)

source/module_esolver/esolver_fp.h

@@ -3,6 +3,10 @@
 
 #include "esolver.h"
 
+#ifndef __MPI
+#include <chrono>
+#endif
+
 //! plane wave basis
 #include "module_basis/module_pw/pw_basis.h"
 
@@ -98,6 +102,12 @@ class ESolver_FP: public ESolver
 
     int pw_rho_flag  = false; ///< flag for pw_rho, 0: not initialized, 1: initialized
 
+    //! the start time of scf iteration
+    #ifdef __MPI
+        double iter_time;
+    #else
+        std::chrono::system_clock::time_point iter_time;
+    #endif
 };
 } // namespace ModuleESolver
 

source/module_esolver/esolver_ks.h

@@ -1,12 +1,6 @@
 #ifndef ESOLVER_KS_H
 #define ESOLVER_KS_H
 
-#ifdef __MPI
-#include <mpi.h>
-#else
-#include <chrono>
-#endif
-
 #include <cstring>
 //#include <fstream>
 
@@ -79,13 +73,6 @@ class ESolver_KS : public ESolver_FP
     //! Electronic wavefunctions
     psi::Psi<T>* psi = nullptr;
 
-    //! the start time of scf iteration
-#ifdef __MPI
-    double iter_time;
-#else
-    std::chrono::system_clock::time_point iter_time;
-#endif
-
     std::string basisname;      //! esolver_ks_lcao.cpp
     double esolver_KS_ne = 0.0; //! number of electrons
     double diag_ethr;           //! the threshold for diagonalization

source/module_esolver/esolver_ks_pw.cpp

@@ -645,14 +645,24 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     //----------------------------------------------------------
     // 3) Print out electronic wavefunctions in pw basis
     //----------------------------------------------------------
-	if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
-	{
-		ModuleIO::write_wfc_pw(GlobalV::KPAR, GlobalV::MY_POOL, GlobalV::MY_RANK, 
-				PARAM.inp.nbands, PARAM.inp.nspin, PARAM.globalv.npol,
-				GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, 
-				PARAM.inp.out_wfc_pw, PARAM.inp.ecutwfc, PARAM.globalv.global_out_dir,
-				this->psi[0], this->kv, this->pw_wfc, GlobalV::ofs_running);
-	}
+    if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
+    {
+        ModuleIO::write_wfc_pw(GlobalV::KPAR,
+                               GlobalV::MY_POOL,
+                               GlobalV::MY_RANK,
+                               PARAM.inp.nbands,
+                               PARAM.inp.nspin,
+                               PARAM.globalv.npol,
+                               GlobalV::RANK_IN_POOL,
+                               GlobalV::NPROC_IN_POOL,
+                               PARAM.inp.out_wfc_pw,
+                               PARAM.inp.ecutwfc,
+                               PARAM.globalv.global_out_dir,
+                               this->psi[0],
+                               this->kv,
+                               this->pw_wfc,
+                               GlobalV::ofs_running);
+    }
 
     //----------------------------------------------------------
     // 4) check if oscillate for delta_spin method
@@ -718,11 +728,6 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const
                               this->pw_rhod->ny,
                               this->pw_rhod->nz,
                               this->pw_rhod->nxyz,
-                              this->kv.get_nks(),
-                              this->kv.isk,
-                              this->kv.wk,
-                              this->pw_big->bz,
-                              this->pw_big->nbz,
                               this->chr.ngmc,
                               &ucell,
                               this->psi,
@@ -731,16 +736,29 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep, const
                               this->ctx,
                               this->Pgrid,
                               PARAM.globalv.global_out_dir,
-                              PARAM.inp.if_separate_k);
+                              PARAM.inp.if_separate_k,
+                              this->kv,
+                              GlobalV::KPAR,
+                              GlobalV::MY_POOL,
+                              &this->chr);
     }
 
-
-	// tmp 2025-05-17, mohan note
-	ModuleIO::write_wfc_pw(GlobalV::KPAR, GlobalV::MY_POOL, GlobalV::MY_RANK, 
-			PARAM.inp.nbands, PARAM.inp.nspin, PARAM.globalv.npol,
-			GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, 
-			PARAM.inp.out_wfc_pw, PARAM.inp.ecutwfc, PARAM.globalv.global_out_dir,
-			this->psi[0], this->kv, this->pw_wfc, GlobalV::ofs_running);
+    // tmp 2025-05-17, mohan note
+    ModuleIO::write_wfc_pw(GlobalV::KPAR,
+                           GlobalV::MY_POOL,
+                           GlobalV::MY_RANK,
+                           PARAM.inp.nbands,
+                           PARAM.inp.nspin,
+                           PARAM.globalv.npol,
+                           GlobalV::RANK_IN_POOL,
+                           GlobalV::NPROC_IN_POOL,
+                           PARAM.inp.out_wfc_pw,
+                           PARAM.inp.ecutwfc,
+                           PARAM.globalv.global_out_dir,
+                           this->psi[0],
+                           this->kv,
+                           this->pw_wfc,
+                           GlobalV::ofs_running);
 
     //------------------------------------------------------------------
     //! 5) calculate Wannier functions in pw basis
@@ -931,7 +949,7 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
     //----------------------------------------------------------
 
     //----------------------------------------------------------
-    //! The write_psi_r_1 interface will be removed in the very 
+    //! The write_psi_r_1 interface will be removed in the very
     //! near future. Don't use it!
     //----------------------------------------------------------
     // if (PARAM.inp.out_wfc_r == 1) // Peize Lin add 2021.11.21
@@ -951,19 +969,15 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
                             this->pw_rhod->ny,
                             this->pw_rhod->nz,
                             this->pw_rhod->nxyz,
-                            this->kv.get_nks(),
-                            this->kv.isk,
-                            this->kv.wk,
-                            this->pw_big->bz,
-                            this->pw_big->nbz,
-                            this->chr.ngmc,
                             &ucell,
                             this->psi,
-                            this->pw_rhod,
                             this->pw_wfc,
                             this->ctx,
                             this->Pgrid,
-                            PARAM.globalv.global_out_dir);
+                            PARAM.globalv.global_out_dir,
+                            this->kv,
+                            GlobalV::KPAR,
+                            GlobalV::MY_POOL);
     }
 
     //----------------------------------------------------------

source/module_esolver/esolver_of.cpp

@@ -164,8 +164,12 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep)
     if (PARAM.inp.of_ml_local_test) this->ml_->localTest(this->chr.rho, this->pw_rho);
 #endif
 
-
     bool conv_esolver = false; // this conv_esolver is added by mohan 20250302 
+#ifdef __MPI
+    this->iter_time = MPI_Wtime();
+#else
+    this->iter_time = std::chrono::system_clock::now();
+#endif
 
     while (true)
     {

source/module_esolver/esolver_of_tool.cpp

@@ -390,38 +390,32 @@ void ESolver_OF::print_info(const bool conv_esolver)
 {
     if (this->iter_ == 0)
     {
-        std::cout << "============================== Running OFDFT "
+        std::cout << " ============================= Running OFDFT "
                      "=============================="
                   << std::endl;
-        std::cout << "Iter        Etot(Ha)          Mu(Ha)      Theta      "
-                     "PotNorm     deltaE(Ha)"
+        std::cout << " ITER       ETOT/eV           EDIFF/eV        EFERMI/eV    POTNORM   TIME/s"
                   << std::endl;
-        // cout << "======================================== Running OFDFT
-        // ========================================" << endl; cout << "Iter
-        // Etot(Ha)          Theta       PotNorm        min/max(den)
-        // min/max(dE/dPhi)" << endl;
     }
-    // ============ used to compare with PROFESS3.0 ================
-    // double minDen = this->chr.rho[0][0];
-    // double maxDen = this->chr.rho[0][0];
-    // double minPot = this->pdEdphi_[0][0];
-    // double maxPot = this->pdEdphi_[0][0];
-    // for (int i = 0; i < this->pw_rho->nrxx; ++i)
-    // {
-    //     if (this->chr.rho[0][i] < minDen) minDen =
-    //     this->chr.rho[0][i]; if (this->chr.rho[0][i] > maxDen)
-    //     maxDen = this->chr.rho[0][i]; if (this->pdEdphi_[0][i] < minPot)
-    //     minPot = this->pdEdphi_[0][i]; if (this->pdEdphi_[0][i] > maxPot)
-    //     maxPot = this->pdEdphi_[0][i];
-    // }
-    std::cout << std::setw(6) << this->iter_ << std::setw(22) << std::scientific << std::setprecision(12)
-              << this->energy_current_ / 2. << std::setw(12) << std::setprecision(3) << this->pelec->eferm.get_efval(0) / 2.
-              << std::setw(12) << this->theta_[0] << std::setw(12) << this->normdLdphi_ << std::setw(12)
-              << (this->energy_current_ - this->energy_last_) / 2. << std::endl;
-    // ============ used to compare with PROFESS3.0 ================
-    // << setw(10) << minDen << "/ " << setw(12) << maxDen
-    // << setw(10) << minPot << "/ " << setw(10) << maxPot << endl;
-    // =============================================================
+
+    std::map<std::string, std::string> prefix_map = {
+        {"cg1", "CG"},
+        {"cg2", "CG"},
+        {"tn", "TN"}
+    };
+    std::string iteration = prefix_map[PARAM.inp.of_method] + std::to_string(this->iter_);
+#ifdef __MPI
+    double duration = (double)(MPI_Wtime() - this->iter_time);
+#else
+    double duration
+        = (std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now() - this->iter_time)).count()
+          / static_cast<double>(1e6);
+#endif
+    std::cout << " " << std::setw(8) << iteration
+              << std::setw(18) << std::scientific << std::setprecision(8) << this->energy_current_ * ModuleBase::Ry_to_eV
+              << std::setw(18) << (this->energy_current_ - this->energy_last_) * ModuleBase::Ry_to_eV
+              << std::setw(13) << std::setprecision(4) << this->pelec->eferm.get_efval(0) * ModuleBase::Ry_to_eV
+              << std::setw(13) << std::setprecision(4) << this->normdLdphi_
+              << std::setw(6) << std::fixed << std::setprecision(2) << duration << std::endl;
 
     GlobalV::ofs_running << std::setprecision(12);
     GlobalV::ofs_running << std::setiosflags(std::ios::right);
@@ -533,5 +527,12 @@ void ESolver_OF::print_info(const bool conv_esolver)
                    /*formats=*/{"%20s", "%20.12f", "%20.12f"}, 0);
     table << titles << energies_Ry << energies_eV;
     GlobalV::ofs_running << table.str() << std::endl;
+
+    // reset the iter_time for the next iteration
+#ifdef __MPI
+    this->iter_time = MPI_Wtime();
+#else
+    this->iter_time = std::chrono::system_clock::now();
+#endif
 }
 } // namespace ModuleESolver