Merge branch 'develop' into Refactor

Author: ErjieWu

docs/community/faq.md

@@ -109,7 +109,7 @@ write(cs_stru, cs_atoms, format='abacus', pp=pp, basis=basis)
 
 ABACUS applies the density difference between two SCF steps (labeled as `DRHO` in the screen output) as the convergence criterion, which is considered as a more robust choice compared with the energy difference. `DRHO` is calculated via `DRHO = |rho(G)-rho_previous(G)|^2`. Note that the energy difference between two SCF steps (labed as `EDIFF`) is also printed out in the screen output.
 
-**4. Why EDIFF is much slower than DRHO?
+**4. Why EDIFF is much slower than DRHO?**
 
 For metaGGA calculations, it is normal because in addition to charge density, kinetic density also needs to be considered in metaGGA calculations. In this case, you can try set `mixing_tau = true`. If you find EDIFF is much slower than DRHO for non-metaGGA calculations, please start a new issue to us.
 

docs/quick_start/easy_install.md

@@ -192,6 +192,8 @@ OMP_NUM_THREADS=4 mpirun -n 4 abacus
 
 In this case, the total thread count is 16.
 
+> Notice: If the MPI library you are using is OpenMPI, which is commonly the case, when you set the number of processes to 1 or 2, OpenMPI will default to `--bind-to core`. This means that no matter how many threads you set, these threads will be restricted to run on 1 or 2 CPU cores. Therefore, setting a higher number of OpenMP threads might result in slower program execution. Hence, when using `mpirun -n` set to 1 or 2, it is recommended to set `--bind-to none` to avoid performance degradation. For example:`OMP_NUM_THREADS=6 mpirun --bind-to none -n 1 abacus`. The detailed binding strategy of OpenMPI can be referred to at https://docs.open-mpi.org/en/v5.0.x/man-openmpi/man1/mpirun.1.html#quick-summary.
+
 ABACUS will try to determine the number of threads used by each process if `OMP_NUM_THREADS` is not set. However, it is **required** to set `OMP_NUM_THREADS` before running `mpirun` to avoid potential performance issues.
 
 Please refer to [hands-on guide](./hands_on.md) for more instructions.

source/module_base/complexarray.cpp

@@ -262,18 +262,4 @@ void point_mult(ComplexArray &in1, ComplexArray &in2, ComplexArray &out){
 		              in1.ptr[i].real() * in2.ptr[i].imag() +
 		              in1.ptr[i].imag() * in2.ptr[i].real());}
 }
-const std::complex <double> &ComplexArray::operator()(const int ind1, const int ind2, const int ind3, const int ind4) const{
-	assert(ind1>=0);	assert(ind1<bound1);
-	assert(ind2>=0);	assert(ind2<bound2);
-	assert(ind3>=0);	assert(ind3<bound3);
-	assert(ind4>=0);	assert(ind4<bound4);
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];}
-std::complex<double>& ComplexArray::operator()(const int ind1,const int ind2,const int ind3,const int ind4){
-	assert(ind1>=0);	assert(ind1<bound1);
-	assert(ind2>=0);	assert(ind2<bound2);
-	assert(ind3>=0);	assert(ind3<bound3);
-	assert(ind4>=0);	assert(ind4<bound4);
-	const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
-	return ptr[ind];}
 }

source/module_base/complexarray.h

@@ -5,6 +5,7 @@
 #include <iostream>
 #include <fstream>
 #include <iomanip>
+#include <cassert>
 
 namespace ModuleBase
 {
@@ -58,11 +59,27 @@ class ComplexArray
 
 	/// overloaded subscript operator for non-const std::complex Array const reference return creates an lvakue
 	std::complex <double> &operator()
-		(const int ind1=0, const int ind2=0, const int ind3=0, const int ind4=0);
+		(const int ind1=0, const int ind2=0, const int ind3=0, const int ind4=0)
+		{
+			assert(ind1>=0);	assert(ind1<bound1);
+			assert(ind2>=0);	assert(ind2<bound2);
+			assert(ind3>=0);	assert(ind3<bound3);
+			assert(ind4>=0);	assert(ind4<bound4);
+			const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+			return ptr[ind];
+		};
 	//  std::complex < double> &operator()(int, int, int, int, int);
 	/// overloaded subscript operator for const std::complex Array const reference return creates an cvakue
 	const std::complex <double> &operator()
-		(const int ind1=0, const int ind2=0, const int ind3=0, const int ind4=0)const;
+		(const int ind1=0, const int ind2=0, const int ind3=0, const int ind4=0) const
+		{
+			assert(ind1>=0);	assert(ind1<bound1);
+			assert(ind2>=0);	assert(ind2<bound2);
+			assert(ind3>=0);	assert(ind3<bound3);
+			assert(ind4>=0);	assert(ind4<bound4);
+			const int ind = ((ind1 * bound2 + ind2) * bound3 + ind3) * bound4 + ind4;
+			return ptr[ind];
+		};
 	//  const std::complex < double> &operator()(int, int, int, int, int)const;
 
 	/****************************************************

source/module_esolver/esolver_fp.cpp

@@ -23,15 +23,16 @@ namespace ModuleESolver
 
 ESolver_FP::ESolver_FP()
 {
-    // pw_rho = new ModuleBase::PW_Basis();
-    // LCAO basis doesn't support GPU acceleration on FFT currently
     std::string fft_device = PARAM.inp.device;
+
+    // LCAO basis doesn't support GPU acceleration on FFT currently
     if(PARAM.inp.basis_type == "lcao")
     {
         fft_device = "cpu";
     }
+
     pw_rho = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
-    if ( PARAM.globalv.double_grid)
+    if (PARAM.globalv.double_grid)
     {
         pw_rhod = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
     }
@@ -44,6 +45,7 @@ ESolver_FP::ESolver_FP()
     pw_big = static_cast<ModulePW::PW_Basis_Big*>(pw_rhod);
     pw_big->setbxyz(PARAM.inp.bx, PARAM.inp.by, PARAM.inp.bz);
     sf.set(pw_rhod, PARAM.inp.nbspline);
+
 }
 
 ESolver_FP::~ESolver_FP()
@@ -140,7 +142,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     // 2) write fermi energy
     ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);
 
-    // 3) update delta rho for charge extrapolation
+    // 3) update delta_rho for charge extrapolation
     CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
 
     if (istep % PARAM.inp.out_interval == 0)
@@ -153,13 +155,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 double* data = nullptr;
                 if (PARAM.inp.dm_to_rho)
                 {
-                    data = this->pelec->charge->rho[is];
-                    this->pw_rhod->real2recip(this->pelec->charge->rho[is], this->pelec->charge->rhog[is]);
+                    data = this->chr.rho[is];
+                    this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]);
                 }
                 else
                 {
-                    data = this->pelec->charge->rho_save[is];
-                    this->pw_rhod->real2recip(this->pelec->charge->rho_save[is], this->pelec->charge->rhog_save[is]);
+                    data = this->chr.rho_save[is];
+                    this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
                 }
                 std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
                 ModuleIO::write_vdata_palgrid(Pgrid,
@@ -176,7 +178,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 {
                     fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube";
                     ModuleIO::write_vdata_palgrid(Pgrid,
-                                                  this->pelec->charge->kin_r_save[is],
+                                                  this->chr.kin_r_save[is],
                                                   is,
                                                   PARAM.inp.nspin,
                                                   istep,
@@ -217,7 +219,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 fn,
                 istep,
                 this->pw_rhod,
-                this->pelec->charge,
+                &this->chr,
                 &(ucell),
                 this->pelec->pot->get_fixed_v(),
                 this->solvent);
@@ -226,11 +228,11 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
         // 6) write ELF
         if (PARAM.inp.out_elf[0] > 0)
         {
-            this->pelec->charge->cal_elf = true;
+            this->chr.cal_elf = true;
             Symmetry_rho srho;
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
-                srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
+                srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
             }
 
             std::string out_dir =PARAM.globalv.global_out_dir;
@@ -242,8 +244,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 out_dir,
                 istep,
                 PARAM.inp.nspin,
-                this->pelec->charge->rho,
-                this->pelec->charge->kin_r,
+                this->chr.rho,
+                this->chr.kin_r,
                 this->pw_rhod,
                 this->Pgrid,
                 &(ucell),
@@ -260,9 +262,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     {
         // only G-vector and K-vector are changed due to the change of lattice
         // vector FFT grids do not change!!
-        pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
-        pw_rho->collect_local_pw();
-        pw_rho->collect_uniqgg();
+        this->pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
+        this->pw_rho->collect_local_pw();
+        this->pw_rho->collect_uniqgg();
 
         if (PARAM.globalv.double_grid)
         {
@@ -292,7 +294,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->CE.update_all_dis(ucell);
         this->CE.extrapolate_charge(&(this->Pgrid),
                                     ucell,
-                                    this->pelec->charge,
+                                    &this->chr,
                                     &(this->sf),
                                     GlobalV::ofs_running,
                                     GlobalV::ofs_warning);
@@ -327,7 +329,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
             std::stringstream ss;
             ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
             ModuleIO::write_vdata_palgrid(this->Pgrid,
-                                          this->pelec->charge->rho[is],
+                                          this->chr.rho[is],
                                           is,
                                           PARAM.inp.nspin,
                                           istep,
@@ -368,8 +370,8 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
         if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
         {
             std::complex<double>** rhog_tot
-                = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rhog : this->pelec->charge->rhog_save;
-            double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rho : this->pelec->charge->rho_save;
+                = (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save;
+            double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save;
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
                 this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]);
@@ -386,12 +388,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
 
             if (XC_Functional::get_ked_flag())
             {
-                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->pelec->charge->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->chr.ngmc, {0.0, 0.0});
                 std::vector<std::complex<double>*> kin_g;
                 for (int is = 0; is < PARAM.inp.nspin; is++)
                 {
-                    kin_g.push_back(kin_g_space.data() + is * this->pelec->charge->ngmc);
-                    this->pw_rhod->real2recip(this->pelec->charge->kin_r_save[is], kin_g[is]);
+                    kin_g.push_back(kin_g_space.data() + is * this->chr.ngmc);
+                    this->pw_rhod->real2recip(this->chr.kin_r_save[is], kin_g[is]);
                 }
                 ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-TAU-DENSITY.restart",
                                      PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local,

source/module_esolver/esolver_fp.h

@@ -2,16 +2,31 @@
 #define ESOLVER_FP_H
 
 #include "esolver.h"
+
+//! plane wave basis
 #include "module_basis/module_pw/pw_basis.h"
+
+//! symmetry analysis
 #include "module_cell/module_symmetry/symmetry.h"
+
+//! electronic states
 #include "module_elecstate/elecstate.h"
+
+//! charge extrapolation
 #include "module_elecstate/module_charge/charge_extra.h"
+
+//! solvation model
 #include "module_hamilt_general/module_surchem/surchem.h"
+
+//! local pseudopotential
 #include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"
+
+//! structure factor related to plane wave basis
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 #include <fstream>
 
+
 //! The First-Principles (FP) Energy Solver Class
 /**
  * This class represents components that needed in
@@ -22,7 +37,7 @@
 
 namespace ModuleESolver
 {
-class ESolver_FP : public ESolver
+class ESolver_FP: public ESolver
 {
   public:
     //! Constructor
@@ -49,39 +64,34 @@ class ESolver_FP : public ESolver
     //! ------------------------------------------------------------------------------
     elecstate::ElecState* pelec = nullptr; ///< Electronic states
 
-    //! ------------------------------------------------------------------------------
+    //! K points in Brillouin zone
+    K_Vectors kv;
 
     //! Electorn charge density
     Charge chr;
 
-    //! Structure factors that used with plane-wave basis set
-    Structure_Factor sf;
-
-    //! K points in Brillouin zone
-    K_Vectors kv;
-
-    //! Plane-wave basis set for charge density
+    //! pw_rho: Plane-wave basis set for charge density
+    //! pw_rhod: same as pw_rho for NCPP. Here 'd' stands for 'dense',
+    //!          dense grid for for uspp, used for ultrasoft augmented charge density.
+    //!          charge density and potential are defined on dense grids,
+    //!          but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
     ModulePW::PW_Basis* pw_rho;
+    ModulePW::PW_Basis* pw_rhod;    //! dense grid for USPP
+    ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class
 
     //! parallel for rho grid
     Parallel_Grid Pgrid;
 
-    //! pointer to local pseudopotential
-    pseudopot_cell_vl locpp;
+    //! Structure factors that used with plane-wave basis set
+    Structure_Factor sf;
 
-    /**
-     * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense'
-     * dense grid for for uspp, used for ultrasoft augmented charge density.
-     * charge density and potential are defined on dense grids,
-     * but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
-     */
-    ModulePW::PW_Basis* pw_rhod;
-    ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class
+    //! local pseudopotentials
+    pseudopot_cell_vl locpp;
 
-    //! Charge extrapolation
+    //! charge extrapolation method
     Charge_Extra CE;
 
-    // solvent model
+    //! solvent model
     surchem solvent;
 };
 } // namespace ModuleESolver