From 193623c80d340b89ef98d97989e00aa72ed705fa Mon Sep 17 00:00:00 2001 From: dyzheng Date: Mon, 30 Sep 2024 16:37:14 +0800 Subject: [PATCH 1/6] Feature: update new version of dav_subspace --- python/pyabacus/src/py_diago_dav_subspace.hpp | 4 +- python/pyabacus/src/py_hsolver.cpp | 6 +- .../pyabacus/src/pyabacus/hsolver/_hsolver.py | 15 +- source/module_elecstate/elecstate.cpp | 25 - source/module_elecstate/elecstate.h | 8 - .../module_elecstate/potentials/pot_local.cpp | 4 + .../module_elecstate/potentials/pot_local.h | 9 + source/module_esolver/esolver_ks_pw.cpp | 10 - source/module_esolver/pw_fun.cpp | 10 - source/module_hsolver/diago_dav_subspace.cpp | 444 +++++------------- source/module_hsolver/diago_dav_subspace.h | 15 +- source/module_hsolver/hsolver_pw.cpp | 47 +- source/module_hsolver/hsolver_pw.h | 6 +- .../module_hsolver/test/test_hsolver_pw.cpp | 4 - source/module_lr/hsolver_lrtd.cpp | 3 +- 15 files changed, 184 insertions(+), 426 deletions(-) diff --git a/python/pyabacus/src/py_diago_dav_subspace.hpp b/python/pyabacus/src/py_diago_dav_subspace.hpp index bd8bbb3e41..2f064787db 100644 --- a/python/pyabacus/src/py_diago_dav_subspace.hpp +++ b/python/pyabacus/src/py_diago_dav_subspace.hpp @@ -106,7 +106,7 @@ class PyDiagoDavSubspace double tol, int max_iter, bool need_subspace, - std::vector is_occupied, + std::vector diag_ethr, bool scf_type, hsolver::diag_comm_info comm_info ) { @@ -143,7 +143,7 @@ class PyDiagoDavSubspace comm_info ); - return obj->diag(hpsi_func, psi, nbasis, eigenvalue, is_occupied, scf_type); + return obj->diag(hpsi_func, psi, nbasis, eigenvalue, diag_ethr.data(), scf_type); } private: diff --git a/python/pyabacus/src/py_hsolver.cpp b/python/pyabacus/src/py_hsolver.cpp index df8f678256..e1c65d876e 100644 --- a/python/pyabacus/src/py_hsolver.cpp +++ b/python/pyabacus/src/py_hsolver.cpp @@ -59,8 +59,8 @@ void bind_hsolver(py::module& m) The maximum number of iterations. need_subspace : bool Whether to use the subspace function. - is_occupied : list[bool] - A list of boolean values indicating whether the band is occupied, + diag_ethr : list[float] + A list of float values indicating the thresholds of each band for the diagonalization, meaning that the corresponding eigenvalue is to be calculated. scf_type : bool Whether to use the SCF type, which is used to determine the @@ -76,7 +76,7 @@ void bind_hsolver(py::module& m) "tol"_a, "max_iter"_a, "need_subspace"_a, - "is_occupied"_a, + "diag_ethr"_a, "scf_type"_a, "comm_info"_a) .def("set_psi", &py_hsolver::PyDiagoDavSubspace::set_psi, R"pbdoc( diff --git a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py index 3dea18ef36..87ad97664e 100644 --- a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py +++ b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py @@ -25,7 +25,7 @@ def dav_subspace( tol: float = 1e-2, max_iter: int = 1000, need_subspace: bool = False, - is_occupied: Union[List[bool], None] = None, + diag_ethr: Union[List[float], None] = None, scf_type: bool = False ) -> Tuple[NDArray[np.float64], NDArray[np.complex128]]: """ A function to diagonalize a matrix using the Davidson-Subspace method. @@ -52,10 +52,8 @@ def dav_subspace( The maximum number of iterations, by default 1000. need_subspace : bool, optional Whether to use subspace function, by default False. - is_occupied : List[bool] | None, optional - The list of occupied bands, by default None. This indicates how many eigenvalues - need to be calculated, starting from the smallest eigenvalue. Only the energy levels - occupied by electrons (occupied) need to be calculated. + diag_ethr : List[float] | None, optional + The list of thresholds of bands, by default None. scf_type : bool, optional Indicates whether the calculation is a self-consistent field (SCF) calculation. If True, the initial precision of eigenvalue calculation can be coarse. @@ -72,8 +70,8 @@ def dav_subspace( if not callable(mvv_op): raise TypeError("mvv_op must be a callable object.") - if is_occupied is None: - is_occupied = [True] * num_eigs + if diag_ethr is None: + diag_ethr = [tol] * num_eigs if init_v.ndim != 1 or init_v.dtype != np.complex128: init_v = init_v.flatten().astype(np.complex128, order='C') @@ -93,7 +91,7 @@ def dav_subspace( tol, max_iter, need_subspace, - is_occupied, + diag_ethr, scf_type, comm_info ) @@ -113,7 +111,6 @@ def davidson( tol: float = 1e-2, max_iter: int = 1000, use_paw: bool = False, - # is_occupied: Union[List[bool], None] = None, # scf_type: bool = False ) -> Tuple[NDArray[np.float64], NDArray[np.complex128]]: """ A function to diagonalize a matrix using the Davidson-Subspace method. diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index 6d77e815a6..26b1623ef6 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -253,31 +253,6 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge this->wg.create(nk_in, GlobalV::NBANDS); } -void set_is_occupied(std::vector& is_occupied, - elecstate::ElecState* pes, - const int i_scf, - const int nk, - const int nband, - const bool diago_full_acc) -{ - if (i_scf != 0 && diago_full_acc == false) - { - for (int i = 0; i < nk; i++) - { - if (pes->klist->wk[i] > 0.0) - { - for (int j = 0; j < nband; j++) - { - if (pes->wg(i, j) / pes->klist->wk[i] < 0.01) - { - is_occupied[i * nband + j] = false; - } - } - } - } - } -}; - } // namespace elecstate diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index 4f1d74683d..624ef355f5 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -194,13 +194,5 @@ class ElecState bool skip_weights = false; }; -// This is an independent function under the elecstate namespace and does not depend on any class. -void set_is_occupied(std::vector& is_occupied, - elecstate::ElecState* pes, - const int i_scf, - const int nk, - const int nband, - const bool diago_full_acc); - } // namespace elecstate #endif diff --git a/source/module_elecstate/potentials/pot_local.cpp b/source/module_elecstate/potentials/pot_local.cpp index 22061b0862..3f9163c132 100644 --- a/source/module_elecstate/potentials/pot_local.cpp +++ b/source/module_elecstate/potentials/pot_local.cpp @@ -8,6 +8,8 @@ namespace elecstate { +double PotLocal::vl_of_0 = 0.0; + //========================================================== // This routine computes the local potential in real space //========================================================== @@ -29,6 +31,8 @@ void PotLocal::cal_fixed_v(double *vl_pseudo // store the local pseudopotential } } + PotLocal::vl_of_0 = vg[0].real(); + // recip2real should be a const function, but now it isn't // a dangerous usage appears here, which should be fix in the future. const_cast(this->rho_basis_)->recip2real(vg, vl_pseudo); diff --git a/source/module_elecstate/potentials/pot_local.h b/source/module_elecstate/potentials/pot_local.h index 327946bc85..9f68dd7946 100644 --- a/source/module_elecstate/potentials/pot_local.h +++ b/source/module_elecstate/potentials/pot_local.h @@ -24,6 +24,15 @@ class PotLocal : public PotBase void cal_fixed_v(double* vl_pseudo) override; + /// @brief get the value of vloc at G=0; + /// @return vl(0) + static double get_vl_of_0() { return vl_of_0; } + + private: + + /// @brief save the value of vloc at G=0; this is a static member because there is only one vl(0) for all instances + static double vl_of_0; + // std::vector vltot; const ModuleBase::matrix* vloc_ = nullptr; // local pseduopotentials diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index af75d8175f..c736d8e239 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -352,15 +352,6 @@ void ESolver_KS_PW::hamilt2density(const int istep, const int iter, c hsolver::DiagoIterAssist::SCF_ITER = iter; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; hsolver::DiagoIterAssist::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax; - - std::vector is_occupied(this->kspw_psi->get_nk() * this->kspw_psi->get_nbands(), true); - - elecstate::set_is_occupied(is_occupied, - this->pelec, - hsolver::DiagoIterAssist::SCF_ITER, - this->kspw_psi->get_nk(), - this->kspw_psi->get_nbands(), - PARAM.inp.diago_full_acc); hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, &this->wf, @@ -383,7 +374,6 @@ void ESolver_KS_PW::hamilt2density(const int istep, const int iter, c this->kspw_psi[0], this->pelec, this->pelec->ekb.c, - is_occupied, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, false); diff --git a/source/module_esolver/pw_fun.cpp b/source/module_esolver/pw_fun.cpp index dd22724598..75dc627795 100644 --- a/source/module_esolver/pw_fun.cpp +++ b/source/module_esolver/pw_fun.cpp @@ -71,15 +71,6 @@ void ESolver_KS_PW::hamilt2estates(const double ethr) hsolver::DiagoIterAssist::need_subspace = false; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - std::vector is_occupied(this->kspw_psi->get_nk() * this->kspw_psi->get_nbands(), true); - - elecstate::set_is_occupied(is_occupied, - this->pelec, - hsolver::DiagoIterAssist::SCF_ITER, - this->kspw_psi->get_nk(), - this->kspw_psi->get_nbands(), - PARAM.inp.diago_full_acc); - hsolver::HSolverPW hsolver_pw_obj(this->pw_wfc, &this->wf, @@ -101,7 +92,6 @@ void ESolver_KS_PW::hamilt2estates(const double ethr) this->kspw_psi[0], this->pelec, this->pelec->ekb.c, - is_occupied, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, true); diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp index bcb9780271..3ecd76ff41 100644 --- a/source/module_hsolver/diago_dav_subspace.cpp +++ b/source/module_hsolver/diago_dav_subspace.cpp @@ -20,8 +20,8 @@ Diago_DavSubspace::Diago_DavSubspace(const std::vector& precond const int& diag_nmax_in, const bool& need_subspace_in, const diag_comm_info& diag_comm_in) - : n_band(nband_in), dim(nbasis_in), nbase_x(nband_in * david_ndim_in), - diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in), precondition(precondition_in) + : precondition(precondition_in), n_band(nband_in), dim(nbasis_in), nbase_x(nband_in * david_ndim_in), + diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in) { this->device = base_device::get_device_type(this->ctx); @@ -59,7 +59,7 @@ Diago_DavSubspace::Diago_DavSubspace(const std::vector& precond if (this->device == base_device::GpuDevice) { resmem_real_op()(this->ctx, this->d_precondition, nbasis_in); - syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition.data(), nbasis_in); + // syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition.data(), nbasis_in); } #endif } @@ -87,7 +87,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, T* psi_in, const int psi_in_dmax, Real* eigenvalue_in_hsolver, - const std::vector& is_occupied) + const double* ethr_band) { ModuleBase::timer::tick("Diago_DavSubspace", "diag_once"); @@ -97,7 +97,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, // convflag[m] = true if the m th band is convergent std::vector convflag(this->n_band, false); - // unconv[m] store the index of the m th unconvergent band + // unconv[m] store the number of the m th unconvergent band std::vector unconv(this->n_band); // the dimension of the reduced psi set @@ -108,9 +108,6 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, ModuleBase::timer::tick("Diago_DavSubspace", "first"); - // copy psi_in to psi_in_iter for the first n_band bands - // NOTE: not the first nbase_x (nband * david_ndim) bands - // since psi_in_iter is zero-initialized, the rest bands are zero for (int m = 0; m < this->n_band; m++) { unconv[m] = m; @@ -122,24 +119,11 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, this->dim); } - // compute h*psi_in_iter - // NOTE: bands after the first n_band should yield zero - hpsi_func(this->psi_in_iter, this->hphi, this->nbase_x, this->dim, 0, this->nbase_x - 1); + hpsi_func(this->psi_in_iter, this->hphi, this->notconv, this->dim, 0, this->notconv - 1); - // at this stage, notconv = n_band and nbase = 0 - // note that nbase of cal_elem is an inout parameter: nbase := nbase + notconv this->cal_elem(this->dim, nbase, this->notconv, this->psi_in_iter, this->hphi, this->hcc, this->scc); - // generalized eigenvalue problem (hcc, scc) for the first n_band bands - this->diag_zhegvx(nbase, - this->n_band, - this->hcc, - this->scc, - this->nbase_x, - &eigenvalue_iter, - this->vcc, - true, - this->is_subspace); + this->diag_zhegvx(nbase, this->notconv, this->hcc, this->scc, this->nbase_x, &eigenvalue_iter, this->vcc); for (size_t m = 0; m < this->n_band; m++) { @@ -154,7 +138,6 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, { dav_iter++; - // nbase is the effective subspace size this->cal_grad(hpsi_func, this->dim, nbase, @@ -167,15 +150,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, this->cal_elem(this->dim, nbase, this->notconv, this->psi_in_iter, this->hphi, this->hcc, this->scc); - this->diag_zhegvx(nbase, - this->n_band, - this->hcc, - this->scc, - this->nbase_x, - &eigenvalue_iter, - this->vcc, - false, - false); + this->diag_zhegvx(nbase, this->n_band, this->hcc, this->scc, this->nbase_x, &eigenvalue_iter, this->vcc); // check convergence and update eigenvalues ModuleBase::timer::tick("Diago_DavSubspace", "check_update"); @@ -183,15 +158,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, this->notconv = 0; for (int m = 0; m < this->n_band; m++) { - if (is_occupied[m]) // always true - { - convflag[m] = (std::abs(eigenvalue_iter[m] - eigenvalue_in_hsolver[m]) < this->diag_thr); - } - else - { - const double empty_ethr = std::max(this->diag_thr * 5.0, 1e-5); - convflag[m] = (std::abs(eigenvalue_iter[m] - eigenvalue_in_hsolver[m]) < empty_ethr); - } + convflag[m] = (std::abs(eigenvalue_iter[m] - eigenvalue_in_hsolver[m]) < ethr_band[m]); if (!convflag[m]) { @@ -230,26 +197,6 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, { // overall convergence or last iteration: exit the iteration - // update this->psi_in_iter according to psi_in - for (size_t i = 0; i < this->n_band; i++) - { - syncmem_complex_op()(this->ctx, - this->ctx, - this->psi_in_iter + i * this->dim, - psi_in + i * psi_in_dmax, - this->dim); - } - - this->refresh(this->dim, - this->n_band, - nbase, - eigenvalue_in_hsolver, - this->psi_in_iter, - this->hphi, - this->hcc, - this->scc, - this->vcc); - ModuleBase::timer::tick("Diago_DavSubspace", "last"); break; } @@ -312,8 +259,6 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, } } - // psi_iter[:, nbase:nbase+notconv] - // = psi_iter[:, :nbase] * vcc[:nbase, :notconv] gemm_op()(this->ctx, 'N', 'N', @@ -321,35 +266,39 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, notconv, nbase, this->one, - psi_iter, + hphi, this->dim, vcc, this->nbase_x, this->zero, - psi_iter + nbase * this->dim, + psi_iter + (nbase) * this->dim, this->dim); - // for psi_iter[:, nbase:nbase+notconv], - // each column is multiplied by the corresponding (minus) eigenvalue - // NOTE: eigenvalue_iter[m] correspond to psi_iter[:, nbase+m] (to be verified) + std::vector e_temp_cpu(nbase, 0); + Real* e_temp_hd = e_temp_cpu.data(); + if(this->device == base_device::GpuDevice) + { + e_temp_hd = nullptr; + resmem_real_op()(this->ctx, e_temp_hd, nbase); + } for (int m = 0; m < notconv; m++) { - - std::vector e_temp_cpu(this->dim, ((-1) * (*eigenvalue_iter)[m])); - + e_temp_cpu.assign(nbase, (-1.0 * (*eigenvalue_iter)[m])); + if (this->device == base_device::GpuDevice) + { + syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, e_temp_hd, e_temp_cpu.data(), nbase); + } vector_mul_vector_op()(this->ctx, - this->dim, - psi_iter + (nbase + m) * this->dim, - psi_iter + (nbase + m) * this->dim, - e_temp_cpu.data()); + nbase, + vcc + m * this->nbase_x, + vcc + m * this->nbase_x, + e_temp_hd); + } + if(this->device == base_device::GpuDevice) + { + delmem_real_op()(this->ctx, e_temp_hd); } - // psi_iter[:, nbase:nbase+notconv] += hphi[:, :nbase] * vcc[:nbase, :notconv] - // - // in terms of input, psi_iter should be - // psi_iter[:, nbase:nbase+notconv] := - // hpsi[:, :nbase] * vcc[:nbase, :notconv] - e[]*psi_iter[:, nbase:nbase+notconv]*vcc[:nbase, :notconv] - // TODO two gemm operations can be combined into one gemm_op()(this->ctx, 'N', 'N', @@ -357,51 +306,46 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, notconv, nbase, this->one, - hphi, + psi_iter, this->dim, vcc, this->nbase_x, this->one, - psi_iter + (nbase) * this->dim, + psi_iter + nbase * this->dim, this->dim); - // FIXME - // (QE dav solver -> isolver = 0) - // In QE, for davidson, h_diag is computed by g2kin + v_of_0 - // this is passed to g_psi, where the residual vector is element-wise - // multiplied by (h_diag - e * s_siag) - // - // There is a macro "TEST_NEW_PRECONDITIONING" in QE, which - // takes x = h_diag - e * s_siag and define - // denm = 0.5 * ( 1 + x + sqrt(1 + (x-1)^2) ) - // which then precondition the residual vector by psi := psi ./ denm - // note that this preconditioning recovers x when x is large, and - // approaches 0.5*(1+sqrt(2)) when x is close to 0. - // - // - // In the current version of ABACUS, "precondition" is simply g2kin, - // and the tested preconditioning of QE is always active. - // // "precondition!!!" std::vector pre(this->dim, 0.0); for (int m = 0; m < notconv; m++) { for (size_t i = 0; i < this->dim; i++) { - //double x = this->precondition[i] - (*eigenvalue_iter)[m]; + // pre[i] = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]); double x = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]); pre[i] = 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0))); } - vector_div_vector_op()(this->ctx, - this->dim, - psi_iter + (nbase + m) * this->dim, - psi_iter + (nbase + m) * this->dim, - pre.data()); +#if defined(__CUDA) || defined(__ROCM) + if (this->device == base_device::GpuDevice) + { + syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, pre.data(), this->dim); + vector_div_vector_op()(this->ctx, + this->dim, + psi_iter + (nbase + m) * this->dim, + psi_iter + (nbase + m) * this->dim, + this->d_precondition); + } + else +#endif + { + vector_div_vector_op()(this->ctx, + this->dim, + psi_iter + (nbase + m) * this->dim, + psi_iter + (nbase + m) * this->dim, + pre.data()); + } } // "normalize!!!" in order to improve numerical stability of subspace diagonalization - - // column-wise normalize psi_iter[:, nbase:nbase+notconv] std::vector psi_norm(notconv, 0.0); for (size_t i = 0; i < notconv; i++) { @@ -409,7 +353,7 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, this->dim, psi_iter + (nbase + i) * this->dim, psi_iter + (nbase + i) * this->dim, - false); // FIXME why reduce is false? + true); assert(psi_norm[i] > 0.0); psi_norm[i] = sqrt(psi_norm[i]); @@ -438,9 +382,6 @@ void Diago_DavSubspace::cal_elem(const int& dim, { ModuleBase::timer::tick("Diago_DavSubspace", "cal_elem"); - // hcc[:, nbase:] = psi_iter.H * hpsi_iter[:, nbase:] - // while psi_iter has nbase_x columns, only the first nbase + notconv columns are - // involved in the matrix multiplication gemm_op()(this->ctx, 'C', 'N', @@ -456,7 +397,6 @@ void Diago_DavSubspace::cal_elem(const int& dim, &hcc[nbase * this->nbase_x], this->nbase_x); - // scc[:, :nbase] := psi_iter.H * psi_iter[:, nbase:] gemm_op()(this->ctx, 'C', 'N', @@ -536,39 +476,6 @@ void Diago_DavSubspace::cal_elem(const int& dim, const size_t last_nbase = nbase; // init: last_nbase = 0 nbase = nbase + notconv; - - // NOTE: nbase is an in & out parameter! nbase is updated to nbase + notconv - - for (size_t i = 0; i < nbase; i++) - { - if (i >= last_nbase) - { - // ensure the diagonal elements of hcc and scc are real - // NOTE: set_real_tocomplex convert a real to itself, and a complex to a new complex - // whose real part is the input's real part and imaginary part is zero - hcc[i * this->nbase_x + i] = set_real_tocomplex(hcc[i * this->nbase_x + i]); - scc[i * this->nbase_x + i] = set_real_tocomplex(scc[i * this->nbase_x + i]); - } - for (size_t j = std::max(i + 1, last_nbase); j < nbase; j++) - { - // ensure the hermicity of hcc and scc - hcc[i * this->nbase_x + j] = get_conj(hcc[j * this->nbase_x + i]); - scc[i * this->nbase_x + j] = get_conj(scc[j * this->nbase_x + i]); - } - } - - // make hcc[nbase:, nbase:] and scc[nbase:, nbase:] to be zero - for (size_t i = nbase; i < this->nbase_x; i++) - { - for (size_t j = nbase; j < this->nbase_x; j++) - { - hcc[i * this->nbase_x + j] = cs.zero; - scc[i * this->nbase_x + j] = cs.zero; - hcc[j * this->nbase_x + i] = cs.zero; - scc[j * this->nbase_x + i] = cs.zero; - } - } - ModuleBase::timer::tick("Diago_DavSubspace", "cal_elem"); return; } @@ -580,18 +487,49 @@ void Diago_DavSubspace::diag_zhegvx(const int& nbase, T* scc, const int& nbase_x, std::vector* eigenvalue_iter, - T* vcc, - bool init, - bool is_subspace) + T* vcc) { ModuleBase::timer::tick("Diago_DavSubspace", "diag_zhegvx"); - if (is_subspace == false) + if (this->diag_comm.rank == 0) { - if (this->diag_comm.rank == 0) + assert(nbase_x >= std::max(1, nbase)); + + if (this->device == base_device::GpuDevice) { - assert(nbase_x >= std::max(1, nbase)); +#if defined(__CUDA) || defined(__ROCM) + Real* eigenvalue_gpu = nullptr; + resmem_real_op()(this->ctx, eigenvalue_gpu, this->nbase_x); + + syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, eigenvalue_gpu, (*eigenvalue_iter).data(), this->nbase_x); + + T* hcc_gpu = nullptr; + T* scc_gpu = nullptr; + T* vcc_gpu = nullptr; + base_device::memory::resize_memory_op()(this->ctx, hcc_gpu, nbase * nbase); + base_device::memory::resize_memory_op()(this->ctx, scc_gpu, nbase * nbase); + base_device::memory::resize_memory_op()(this->ctx, vcc_gpu, nbase * nbase); + for(int i=0;i()(this->ctx, this->ctx, hcc_gpu + i * nbase, hcc + i * nbase_x, nbase); + base_device::memory::synchronize_memory_op()(this->ctx, this->ctx, scc_gpu + i * nbase, scc + i * nbase_x, nbase); + } + dngvd_op()(this->ctx, nbase, nbase, hcc_gpu, scc_gpu, eigenvalue_gpu, vcc_gpu); + for(int i=0;i()(this->ctx, this->ctx, vcc + i * nbase_x, vcc_gpu + i * nbase, nbase); + } + delmem_complex_op()(this->ctx, hcc_gpu); + delmem_complex_op()(this->ctx, scc_gpu); + delmem_complex_op()(this->ctx, vcc_gpu); + + syncmem_var_d2h_op()(this->cpu_ctx, this->ctx, (*eigenvalue_iter).data(), eigenvalue_gpu, this->nbase_x); + delmem_real_op()(this->ctx, eigenvalue_gpu); +#endif + } + else + { std::vector> h_diag(nbase, std::vector(nbase, cs.zero)); std::vector> s_diag(nbase, std::vector(nbase, cs.zero)); @@ -603,56 +541,14 @@ void Diago_DavSubspace::diag_zhegvx(const int& nbase, s_diag[i][j] = scc[i * this->nbase_x + j]; } } - - if (this->device == base_device::GpuDevice) - { -#if defined(__CUDA) || defined(__ROCM) - Real* eigenvalue_gpu = nullptr; - resmem_real_op()(this->ctx, eigenvalue_gpu, this->nbase_x); - - syncmem_var_h2d_op()(this->ctx, - this->cpu_ctx, - eigenvalue_gpu, - (*eigenvalue_iter).data(), - this->nbase_x); - - dnevx_op()(this->ctx, nbase, this->nbase_x, this->hcc, nband, eigenvalue_gpu, this->vcc); - - syncmem_var_d2h_op()(this->cpu_ctx, - this->ctx, - (*eigenvalue_iter).data(), - eigenvalue_gpu, - this->nbase_x); - - delmem_real_op()(this->ctx, eigenvalue_gpu); -#endif - } - else - { - if (init) - { - dnevx_op()(this->ctx, - nbase, - this->nbase_x, - this->hcc, - nband, - (*eigenvalue_iter).data(), - this->vcc); - } - else - { - - dngvx_op()(this->ctx, - nbase, - this->nbase_x, - this->hcc, - this->scc, - nband, - (*eigenvalue_iter).data(), - this->vcc); - } - } - + dngvx_op()(this->ctx, + nbase, + this->nbase_x, + this->hcc, + this->scc, + nband, + (*eigenvalue_iter).data(), + this->vcc); // reset: for (size_t i = 0; i < nbase; i++) { @@ -671,35 +567,19 @@ void Diago_DavSubspace::diag_zhegvx(const int& nbase, } } } - -#ifdef __MPI - if (this->diag_comm.nproc > 1) - { - // vcc: nbase * nband - for (int i = 0; i < nband; i++) - { - MPI_Bcast(&vcc[i * this->nbase_x], nbase, MPI_DOUBLE_COMPLEX, 0, this->diag_comm.comm); - } - MPI_Bcast((*eigenvalue_iter).data(), nband, MPI_DOUBLE, 0, this->diag_comm.comm); - } -#endif } - else if (is_subspace == true) - { - for (size_t m = 0; m < nband; m++) - { - (*eigenvalue_iter)[m] = get_real(hcc[m * this->nbase_x + m]); - - vcc[m * this->nbase_x + m] = set_real_tocomplex(1.0); - } #ifdef __MPI - if (this->diag_comm.nproc > 1) + if (this->diag_comm.nproc > 1) + { + // vcc: nbase * nband + for (int i = 0; i < nband; i++) { - MPI_Bcast((*eigenvalue_iter).data(), this->n_band, MPI_DOUBLE, 0, this->diag_comm.comm); + MPI_Bcast(&vcc[i * this->nbase_x], nbase, MPI_DOUBLE_COMPLEX, 0, this->diag_comm.comm); } -#endif + MPI_Bcast((*eigenvalue_iter).data(), nband, MPI_DOUBLE, 0, this->diag_comm.comm); } +#endif ModuleBase::timer::tick("Diago_DavSubspace", "diag_zhegvx"); return; @@ -805,7 +685,7 @@ int Diago_DavSubspace::diag(const HPsiFunc& hpsi_func, T* psi_in, const int psi_in_dmax, Real* eigenvalue_in_hsolver, - const std::vector& is_occupied, + const double* ethr_band, const bool& scf_type) { /// record the times of trying iterative diagonalization @@ -816,13 +696,8 @@ int Diago_DavSubspace::diag(const HPsiFunc& hpsi_func, do { - //printf("enter diag... is_subspace = %d, ntry = %d\n", this->is_subspace, ntry); - if (this->is_subspace || ntry > 0) - { - this->diagH_subspace(psi_in, eigenvalue_in_hsolver, hpsi_func, this->n_band, this->dim, psi_in_dmax); - } - sum_iter += this->diag_once(hpsi_func, psi_in, psi_in_dmax, eigenvalue_in_hsolver, is_occupied); + sum_iter += this->diag_once(hpsi_func, psi_in, psi_in_dmax, eigenvalue_in_hsolver, ethr_band); ++ntry; @@ -849,113 +724,12 @@ bool Diago_DavSubspace::test_exit_cond(const int& ntry, const int& no // In non-self consistent calculation, do until totally converged. const bool f2 = ((!scf && (notconv > 0))); - // if self consistent calculation, if not converged > 5, - // using diagH_subspace and cg method again. ntry++ + // if self consistent calculation, if not converged > 5 const bool f3 = ((scf && (notconv > 5))); return (f1 && (f2 || f3)); } -template -void Diago_DavSubspace::diagH_subspace(T* psi_pointer, // [in] & [out] wavefunction - Real* en, // [out] eigenvalues - const HPsiFunc hpsi_func, - const int n_band, - const int dmin, - const int dmax) -{ - ModuleBase::TITLE("Diago_DavSubspace", "subspace"); - ModuleBase::timer::tick("Diago_DavSubspace", "subspace"); - - const int nstart = n_band; - assert(n_band > 0); - - T* hcc = nullptr; - T* scc = nullptr; - T* vcc = nullptr; - resmem_complex_op()(ctx, hcc, nstart * nstart, "DAV::hcc"); - resmem_complex_op()(ctx, scc, nstart * nstart, "DAV::scc"); - resmem_complex_op()(ctx, vcc, nstart * nstart, "DAV::vcc"); - setmem_complex_op()(ctx, hcc, 0, nstart * nstart); - setmem_complex_op()(ctx, scc, 0, nstart * nstart); - setmem_complex_op()(ctx, vcc, 0, nstart * nstart); - - T* hphi = nullptr; - resmem_complex_op()(ctx, hphi, nstart * dmax, "DAV::hphi"); - setmem_complex_op()(ctx, hphi, 0, nstart * dmax); - - { - // do hPsi for all bands - hpsi_func(psi_pointer, hphi, n_band, dmax, 0, nstart - 1); - - gemm_op()(ctx, - 'C', - 'N', - nstart, - nstart, - dmin, - this->one, - psi_pointer, - dmax, - hphi, - dmax, - this->zero, - hcc, - nstart); - - gemm_op()(ctx, - 'C', - 'N', - nstart, - nstart, - dmin, - this->one, - psi_pointer, - dmax, - psi_pointer, - dmax, - this->zero, - scc, - nstart); - } - - if (GlobalV::NPROC_IN_POOL > 1) - { - Parallel_Reduce::reduce_pool(hcc, nstart * nstart); - Parallel_Reduce::reduce_pool(scc, nstart * nstart); - } - - // after generation of H and S matrix, diag them - DiagoIterAssist::diagH_LAPACK(nstart, n_band, hcc, scc, nstart, en, vcc); - - { // code block to calculate evc - gemm_op()(ctx, - 'N', - 'N', - dmin, - n_band, - nstart, - this->one, - psi_pointer, // dmin * nstart - dmax, - vcc, // nstart * n_band - nstart, - this->zero, - hphi, - dmin); - } - - matrixSetToAnother()(ctx, n_band, hphi, dmin, psi_pointer, dmax); - - delmem_complex_op()(ctx, hphi); - - delmem_complex_op()(ctx, hcc); - delmem_complex_op()(ctx, scc); - delmem_complex_op()(ctx, vcc); - - ModuleBase::timer::tick("Diago_DavSubspace", "diagH_subspace"); -} - namespace hsolver { diff --git a/source/module_hsolver/diago_dav_subspace.h b/source/module_hsolver/diago_dav_subspace.h index 78553f0174..17cdc1263c 100644 --- a/source/module_hsolver/diago_dav_subspace.h +++ b/source/module_hsolver/diago_dav_subspace.h @@ -37,7 +37,7 @@ class Diago_DavSubspace : public DiagH T* psi_in, const int psi_in_dmax, Real* eigenvalue_in, - const std::vector& is_occupied, + const double* ethr_band, const bool& scf_type); private: @@ -110,13 +110,6 @@ class Diago_DavSubspace : public DiagH T* scc, T* vcc); - void diagH_subspace(T* psi_pointer, // [in] & [out] wavefunction - Real* en, // [out] eigenvalues - const HPsiFunc hpsi_func, - const int n_band, - const int dmin, - const int dmax); - // void diagH_LAPACK(const int nstart, // const int nbands, // const T* hcc, @@ -131,15 +124,13 @@ class Diago_DavSubspace : public DiagH T* scc, const int& nbase_x, std::vector* eigenvalue_iter, - T* vcc, - bool init, - bool is_subspace); + T* vcc); int diag_once(const HPsiFunc& hpsi_func, T* psi_in, const int psi_in_dmax, Real* eigenvalue_in, - const std::vector& is_occupied); + const double* ethr_band); bool test_exit_cond(const int& ntry, const int& notconv, const bool& scf); diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 7b8e5fc0a6..14d2453453 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -15,6 +15,7 @@ #include "module_hsolver/diago_dav_subspace.h" #include "module_hsolver/diago_david.h" #include "module_hsolver/diago_iter_assist.h" +#include "module_elecstate/potentials/pot_local.h" #include #include @@ -208,12 +209,43 @@ void HSolverPW::paw_func_after_kloop(psi::Psi& psi, elecst #endif +template +void HSolverPW::cal_ethr_band(const double& wk, const double* wg, const double& ethr, std::vector& ethrs) +{ + if(wk > 0.0) + { + const double ethr_unocc = std::max(1e-5, ethr); + for (int i = 0; i < ethrs.size(); i++) + { + double band_weight = wg[i] / wk; + if (band_weight > 1e-2) + { + ethrs[i] = ethr; + } + else if(band_weight > 1e-5) + {// similar energy difference for difsferent bands + ethrs[i] = std::min(ethr_unocc, ethr / band_weight); + } + else + { + ethrs[i] = ethr_unocc; + } + } + } + else + { + for (int i = 0; i < ethrs.size(); i++) + { + ethrs[i] = ethr; + } + } +} + template void HSolverPW::solve(hamilt::Hamilt* pHamilt, psi::Psi& psi, elecstate::ElecState* pes, double* out_eigenvalues, - const std::vector& is_occupied_in, const int rank_in_pool_in, const int nproc_in_pool_in, const bool skip_charge) @@ -233,8 +265,8 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // prepare for the precondition of diagonalization std::vector precondition(psi.get_nbasis(), 0.0); - // std::vector eigenvalues(pes->ekb.nr * pes->ekb.nc, 0.0); std::vector eigenvalues(this->wfc_basis->nks * psi.get_nbands(), 0.0); + ethr_band.resize(psi.get_nbands(), DiagoIterAssist::PW_DIAG_THR); /// Loop over k points for solve Hamiltonian to charge density for (int ik = 0; ik < this->wfc_basis->nks; ++ik) @@ -250,6 +282,10 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // template add precondition calculating here update_precondition(precondition, ik, this->wfc_basis->npwk[ik]); + this->cal_ethr_band(pes->klist->wk[ik], + &pes->wg(ik, 0), + DiagoIterAssist::PW_DIAG_THR, + ethr_band); #ifdef USE_PAW this->call_paw_cell_set_currentk(ik); @@ -454,7 +490,6 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, ModuleBase::timer::tick("DavSubspace", "hpsi_func"); }; bool scf = this->calculation_type == "nscf" ? false : true; - const std::vector is_occupied(psi.get_nbands(), true); Diago_DavSubspace dav_subspace(pre_condition, psi.get_nbands(), @@ -466,8 +501,8 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, this->need_subspace, comm_info); - DiagoIterAssist::avg_iter += static_cast( - dav_subspace.diag(hpsi_func, psi.get_pointer(), psi.get_nbasis(), eigenvalue, is_occupied, scf)); + DiagoIterAssist::avg_iter + += static_cast(dav_subspace.diag(hpsi_func, psi.get_pointer(), psi.get_nbasis(), eigenvalue, this->ethr_band.data(), scf)); } else if (this->method == "dav") { @@ -574,7 +609,7 @@ void HSolverPW::update_precondition(std::vector& h_diag, const if (this->method == "dav_subspace") { - h_diag[ig] = g2kin; + h_diag[ig] = g2kin + Real(elecstate::PotLocal::get_vl_of_0()); } else { diff --git a/source/module_hsolver/hsolver_pw.h b/source/module_hsolver/hsolver_pw.h index 62bbc5120b..72c06b67dc 100644 --- a/source/module_hsolver/hsolver_pw.h +++ b/source/module_hsolver/hsolver_pw.h @@ -48,7 +48,6 @@ class HSolverPW psi::Psi& psi, elecstate::ElecState* pes, double* out_eigenvalues, - const std::vector& is_occupied_in, const int rank_in_pool_in, const int nproc_in_pool_in, const bool skip_charge); @@ -91,6 +90,11 @@ class HSolverPW int rank_in_pool = 0; int nproc_in_pool = 1; + /// @brief calculate the threshold for iterative-diagonalization for each band + void cal_ethr_band(const double& wk, const double* wg, const double& ethr, std::vector& ethrs); + + std::vector ethr_band; + #ifdef USE_PAW void paw_func_in_kloop(const int ik); diff --git a/source/module_hsolver/test/test_hsolver_pw.cpp b/source/module_hsolver/test/test_hsolver_pw.cpp index 2d8d5b4f07..939393f3d3 100644 --- a/source/module_hsolver/test/test_hsolver_pw.cpp +++ b/source/module_hsolver/test/test_hsolver_pw.cpp @@ -93,13 +93,10 @@ TEST_F(TestHSolverPW, solve) { EXPECT_EQ(this->hs_f.initialed_psi, false); EXPECT_EQ(this->hs_d.initialed_psi, false); - std::vector is_occupied(1 * 2, true); - this->hs_f.solve(&hamilt_test_f, psi_test_cf, &elecstate_test, elecstate_test.ekb.c, - is_occupied, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, @@ -118,7 +115,6 @@ TEST_F(TestHSolverPW, solve) { psi_test_cd, &elecstate_test, elecstate_test.ekb.c, - is_occupied, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, diff --git a/source/module_lr/hsolver_lrtd.cpp b/source/module_lr/hsolver_lrtd.cpp index 0cc95fb8f6..2e0f8cc769 100644 --- a/source/module_lr/hsolver_lrtd.cpp +++ b/source/module_lr/hsolver_lrtd.cpp @@ -148,12 +148,13 @@ namespace LR nband_in); }; + std::vector ethr_band(psi_k1_dav.get_nbands(), this->diag_ethr); hsolver::DiagoIterAssist::avg_iter += static_cast(dav_subspace.diag( hpsi_func, psi_k1_dav.get_pointer(), psi_k1_dav.get_nbasis(), eigenvalue.data(), - std::vector(psi_k1_dav.get_nbands(), true), + ethr_band.data(), false /*scf*/)); } else {throw std::runtime_error("HSolverLR::solve: method not implemented");} From a11d1feae80d189cf7b07a1fdaa811402a5c6b14 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Mon, 30 Sep 2024 17:21:52 +0800 Subject: [PATCH 2/6] Fix: UnitTest error --- source/module_hsolver/test/hsolver_supplementary_mock.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index 95c69e555e..c1e81aeec9 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -1,5 +1,6 @@ #pragma once #include "module_elecstate/elecstate.h" +#include "module_elecstate/potentials/pot_local.h" namespace elecstate { @@ -57,6 +58,8 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge return; } +double PotLocal::vl_of_0 = 0.0; + } // namespace elecstate From d4f570da06415e07219c30c14bb908e4cecd2f68 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Mon, 30 Sep 2024 10:28:36 +0000 Subject: [PATCH 3/6] [pre-commit.ci lite] apply automatic fixes --- source/module_lr/hsolver_lrtd.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/module_lr/hsolver_lrtd.cpp b/source/module_lr/hsolver_lrtd.cpp index 2e0f8cc769..6f6ae86a6f 100644 --- a/source/module_lr/hsolver_lrtd.cpp +++ b/source/module_lr/hsolver_lrtd.cpp @@ -14,7 +14,8 @@ namespace LR inline void print_eigs(const std::vector& eigs, const std::string& label = "", const double factor = 1.0) { std::cout << label << std::endl; - for (auto& e : eigs)std::cout << e * factor << " "; + for (auto& e : eigs) {std::cout << e * factor << " "; +} std::cout << std::endl; } template From 7b8d46ff10d30690fcbe5e850f658cc06c139c35 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 8 Oct 2024 16:07:38 +0800 Subject: [PATCH 4/6] Fix: UnitTest HSolver_PW --- source/module_hsolver/hsolver_pw.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 14d2453453..3c581e77f2 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -214,6 +214,8 @@ void HSolverPW::cal_ethr_band(const double& wk, const double* wg, con { if(wk > 0.0) { + // Note: the idea of threshold for unoccupied bands (1e-5) comes from QE + // In ABACUS, We applied a smoothing process to this truncation to avoid abrupt changes in energy errors between different bands. const double ethr_unocc = std::max(1e-5, ethr); for (int i = 0; i < ethrs.size(); i++) { @@ -223,7 +225,7 @@ void HSolverPW::cal_ethr_band(const double& wk, const double* wg, con ethrs[i] = ethr; } else if(band_weight > 1e-5) - {// similar energy difference for difsferent bands + {// similar energy difference for different bands when band_weight in range [1e-5, 1e-2] ethrs[i] = std::min(ethr_unocc, ethr / band_weight); } else @@ -282,10 +284,16 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, // template add precondition calculating here update_precondition(precondition, ik, this->wfc_basis->npwk[ik]); - this->cal_ethr_band(pes->klist->wk[ik], + + // only dav_subspace method used smooth threshold for all bands now, + // for other methods, this trick can be added in the future to accelerate calculation without accuracy loss. + if (this->method == "dav_subspace") + { + this->cal_ethr_band(pes->klist->wk[ik], &pes->wg(ik, 0), DiagoIterAssist::PW_DIAG_THR, ethr_band); + } #ifdef USE_PAW this->call_paw_cell_set_currentk(ik); From 7f1aa92ee2b8fa26e5d3b97829085aedd1c9506a Mon Sep 17 00:00:00 2001 From: dyzheng Date: Fri, 11 Oct 2024 16:26:21 +0800 Subject: [PATCH 5/6] Fix: comments in PR --- docs/advanced/input_files/input-main.md | 9 +-------- examples/lr-tddft/lcao_Si2/INPUT | 1 - .../module_elecstate/potentials/pot_local.cpp | 8 +++++--- source/module_elecstate/potentials/pot_local.h | 11 ++++------- .../module_elecstate/potentials/potential_new.h | 10 ++++++++++ .../potentials/potential_types.cpp | 2 +- source/module_hsolver/hsolver_pw.cpp | 17 ++++++++++------- source/module_hsolver/hsolver_pw.h | 2 +- source/module_hsolver/hsolver_pw_sdft.cpp | 2 +- .../test/hsolver_supplementary_mock.h | 3 --- source/module_io/read_input_item_system.cpp | 14 -------------- source/module_parameter/input_parameter.h | 1 - tests/integrate/102_PW_DS_davsubspace/INPUT | 1 - 13 files changed, 33 insertions(+), 48 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index b332f1fc80..fb5deed5f6 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -39,7 +39,6 @@ - [pw\_diag\_thr](#pw_diag_thr) - [pw\_diag\_nmax](#pw_diag_nmax) - [pw\_diag\_ndim](#pw_diag_ndim) - - [diago\_full\_acc](#diago_full_acc) - [erf\_ecut](#erf_ecut) - [fft\_mode](#fft_mode) - [erf\_height](#erf_height) @@ -779,12 +778,6 @@ These variables are used to control the plane wave related parameters. - **Description**: Only useful when you use `ks_solver = dav` or `ks_solver = dav_subspace`. It indicates dimension of workspace(number of wavefunction packets, at least 2 needed) for the Davidson method. A larger value may yield a smaller number of iterations in the algorithm but uses more memory and more CPU time in subspace diagonalization. - **Default**: 4 -### diago_full_acc - -- **Type**: bool -- **Description**: Only useful when you use `ks_solver = dav_subspace`. If `TRUE`, all the empty states are diagonalized at the same level of accuracy of the occupied ones. Otherwise the empty states are diagonalized using a larger threshold (10-5) (this should not affect total energy, forces, and other ground-state properties). -- **Default**: false - ### erf_ecut - **Type**: Real @@ -925,7 +918,7 @@ calculations. - **cg**: cg method. - **bpcg**: bpcg method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. - **dav**: the Davidson algorithm. - - **dav_subspace**: subspace Davidson algorithm + - **dav_subspace**: Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. For atomic orbitals basis, diff --git a/examples/lr-tddft/lcao_Si2/INPUT b/examples/lr-tddft/lcao_Si2/INPUT index cac5e5cca3..19ea04fe1a 100644 --- a/examples/lr-tddft/lcao_Si2/INPUT +++ b/examples/lr-tddft/lcao_Si2/INPUT @@ -37,4 +37,3 @@ out_alllog 1 nvirt 19 abs_wavelen_range 100 175 -#diago_full_acc 1 diff --git a/source/module_elecstate/potentials/pot_local.cpp b/source/module_elecstate/potentials/pot_local.cpp index 3f9163c132..b387fc0328 100644 --- a/source/module_elecstate/potentials/pot_local.cpp +++ b/source/module_elecstate/potentials/pot_local.cpp @@ -8,8 +8,6 @@ namespace elecstate { -double PotLocal::vl_of_0 = 0.0; - //========================================================== // This routine computes the local potential in real space //========================================================== @@ -31,7 +29,11 @@ void PotLocal::cal_fixed_v(double *vl_pseudo // store the local pseudopotential } } - PotLocal::vl_of_0 = vg[0].real(); + /// save the V_local at G=0 + if(this->rho_basis_->npw > 0) + { + *vl_of_0_ = vg[0].real(); + } // recip2real should be a const function, but now it isn't // a dangerous usage appears here, which should be fix in the future. diff --git a/source/module_elecstate/potentials/pot_local.h b/source/module_elecstate/potentials/pot_local.h index 9f68dd7946..df31ac471b 100644 --- a/source/module_elecstate/potentials/pot_local.h +++ b/source/module_elecstate/potentials/pot_local.h @@ -12,8 +12,9 @@ class PotLocal : public PotBase public: PotLocal(const ModuleBase::matrix* vloc_in, // local pseduopotentials const ModuleBase::ComplexMatrix* sf_in, - const ModulePW::PW_Basis* rho_basis_in) - : vloc_(vloc_in), sf_(sf_in) + const ModulePW::PW_Basis* rho_basis_in, + double& vl_of_0) + : vloc_(vloc_in), sf_(sf_in), vl_of_0_(&vl_of_0) { assert(this->vloc_->nr == this->sf_->nr); this->rho_basis_ = rho_basis_in; @@ -24,14 +25,10 @@ class PotLocal : public PotBase void cal_fixed_v(double* vl_pseudo) override; - /// @brief get the value of vloc at G=0; - /// @return vl(0) - static double get_vl_of_0() { return vl_of_0; } - private: /// @brief save the value of vloc at G=0; this is a static member because there is only one vl(0) for all instances - static double vl_of_0; + double* vl_of_0_ = nullptr; // std::vector vltot; diff --git a/source/module_elecstate/potentials/potential_new.h b/source/module_elecstate/potentials/potential_new.h index f3ef3b4841..d1f8e41f78 100644 --- a/source/module_elecstate/potentials/potential_new.h +++ b/source/module_elecstate/potentials/potential_new.h @@ -169,6 +169,14 @@ class Potential : public PotBase return this->v_effective_fixed.data(); } + + /// @brief get the value of vloc at G=0; + /// @return vl(0) + double get_vl_of_0() const + { + return this->vl_of_0; + } + private: void cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff) override; void cal_fixed_v(double* vl_pseudo) override; @@ -196,6 +204,8 @@ class Potential : public PotBase double* etxc_ = nullptr; double* vtxc_ = nullptr; + double vl_of_0 = 0.0; + std::vector components; const UnitCell* ucell_ = nullptr; diff --git a/source/module_elecstate/potentials/potential_types.cpp b/source/module_elecstate/potentials/potential_types.cpp index 36ac62040d..bc15110798 100644 --- a/source/module_elecstate/potentials/potential_types.cpp +++ b/source/module_elecstate/potentials/potential_types.cpp @@ -27,7 +27,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) { if(!PARAM.inp.use_paw) { - return new PotLocal(this->vloc_, &(this->structure_factors_->strucFac), this->rho_basis_); + return new PotLocal(this->vloc_, &(this->structure_factors_->strucFac), this->rho_basis_, this->vl_of_0); } else { diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 66d17a37cb..43137807f3 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -15,7 +15,6 @@ #include "module_hsolver/diago_dav_subspace.h" #include "module_hsolver/diago_david.h" #include "module_hsolver/diago_iter_assist.h" -#include "module_elecstate/potentials/pot_local.h" #include #include @@ -212,19 +211,23 @@ void HSolverPW::paw_func_after_kloop(psi::Psi& psi, elecst template void HSolverPW::cal_ethr_band(const double& wk, const double* wg, const double& ethr, std::vector& ethrs) { + // threshold for classifying occupied and unoccupied bands + const double occ_threshold = 1e-2; + // diagonalization threshold limitation for unoccupied bands + const double ethr_limit = 1e-5; if(wk > 0.0) { // Note: the idea of threshold for unoccupied bands (1e-5) comes from QE // In ABACUS, We applied a smoothing process to this truncation to avoid abrupt changes in energy errors between different bands. - const double ethr_unocc = std::max(1e-5, ethr); + const double ethr_unocc = std::max(ethr_limit, ethr); for (int i = 0; i < ethrs.size(); i++) { double band_weight = wg[i] / wk; - if (band_weight > 1e-2) + if (band_weight > occ_threshold) { ethrs[i] = ethr; } - else if(band_weight > 1e-5) + else if(band_weight > ethr_limit) {// similar energy difference for different bands when band_weight in range [1e-5, 1e-2] ethrs[i] = std::min(ethr_unocc, ethr / band_weight); } @@ -283,7 +286,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, this->updatePsiK(pHamilt, psi, ik); // template add precondition calculating here - update_precondition(precondition, ik, this->wfc_basis->npwk[ik]); + update_precondition(precondition, ik, this->wfc_basis->npwk[ik], Real(pes->pot->get_vl_of_0())); // only dav_subspace method used smooth threshold for all bands now, // for other methods, this trick can be added in the future to accelerate calculation without accuracy loss. @@ -589,7 +592,7 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, } template -void HSolverPW::update_precondition(std::vector& h_diag, const int ik, const int npw) +void HSolverPW::update_precondition(std::vector& h_diag, const int ik, const int npw, const Real vl_of_0) { h_diag.assign(h_diag.size(), 1.0); int precondition_type = 2; @@ -616,7 +619,7 @@ void HSolverPW::update_precondition(std::vector& h_diag, const if (this->method == "dav_subspace") { - h_diag[ig] = g2kin + Real(elecstate::PotLocal::get_vl_of_0()); + h_diag[ig] = g2kin + vl_of_0; } else { diff --git a/source/module_hsolver/hsolver_pw.h b/source/module_hsolver/hsolver_pw.h index 72c06b67dc..a8ca9f5b16 100644 --- a/source/module_hsolver/hsolver_pw.h +++ b/source/module_hsolver/hsolver_pw.h @@ -63,7 +63,7 @@ class HSolverPW void updatePsiK(hamilt::Hamilt* pHamilt, psi::Psi& psi, const int ik); // calculate the precondition array for diagonalization in PW base - void update_precondition(std::vector& h_diag, const int ik, const int npw); + void update_precondition(std::vector& h_diag, const int ik, const int npw, const Real vl_of_0); void output_iterInfo(); diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp index 4f4c858cf4..0318d6c5fb 100644 --- a/source/module_hsolver/hsolver_pw_sdft.cpp +++ b/source/module_hsolver/hsolver_pw_sdft.cpp @@ -44,7 +44,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt>* pHamilt, { this->updatePsiK(pHamilt, psi, ik); // template add precondition calculating here - update_precondition(precondition, ik, this->wfc_basis->npwk[ik]); + update_precondition(precondition, ik, this->wfc_basis->npwk[ik], pes->pot->get_vl_of_0()); /// solve eigenvector and eigenvalue for H(k) double* p_eigenvalues = &(pes->ekb(ik, 0)); this->hamiltSolvePsiK(pHamilt, psi, precondition, p_eigenvalues); diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index c1e81aeec9..95c69e555e 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -1,6 +1,5 @@ #pragma once #include "module_elecstate/elecstate.h" -#include "module_elecstate/potentials/pot_local.h" namespace elecstate { @@ -58,8 +57,6 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge return; } -double PotLocal::vl_of_0 = 0.0; - } // namespace elecstate diff --git a/source/module_io/read_input_item_system.cpp b/source/module_io/read_input_item_system.cpp index 32ff4bccef..8e36273d2b 100644 --- a/source/module_io/read_input_item_system.cpp +++ b/source/module_io/read_input_item_system.cpp @@ -442,20 +442,6 @@ void ReadInput::item_system() read_sync_int(input.fft_mode); this->add_item(item); } - { - Input_Item item("diago_full_acc"); - item.annotation = "all the empty states are diagonalized"; - /** - * @brief diago_full_acc - * If .TRUE. all the empty states are diagonalized at the same level of - * accuracy of the occupied ones. Otherwise the empty states are - * diagonalized using a larger threshold (this should not affect total - * energy, forces, and other ground-state properties). - * - */ - read_sync_bool(input.diago_full_acc); - this->add_item(item); - } { Input_Item item("init_wfc"); item.annotation = "start wave functions are from 'atomic', " diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 57bcab403d..d984b1288c 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -44,7 +44,6 @@ struct Input_para double erf_height = 0; ///< the height of the energy step for reciprocal vectors double erf_sigma = 0.1; ///< the width of the energy step for reciprocal vectors int fft_mode = 0; ///< fftw mode 0: estimate, 1: measure, 2: patient, 3: exhaustive - bool diago_full_acc = false; ///< all the empty states are diagonalized std::string init_wfc = "atomic"; ///< "file","atomic","random" bool psi_initializer = false; ///< whether use psi_initializer to initialize wavefunctions int pw_seed = 0; ///< random seed for initializing wave functions diff --git a/tests/integrate/102_PW_DS_davsubspace/INPUT b/tests/integrate/102_PW_DS_davsubspace/INPUT index 1a834c95b6..18f7dc5714 100644 --- a/tests/integrate/102_PW_DS_davsubspace/INPUT +++ b/tests/integrate/102_PW_DS_davsubspace/INPUT @@ -25,4 +25,3 @@ mixing_type plain mixing_beta 0.5 ks_solver dav_subspace -diago_full_acc 1 From 5ddc59d06f7f3e2c1c5f125f600e8c87e9572388 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sat, 12 Oct 2024 18:22:30 +0800 Subject: [PATCH 6/6] Fix: error in UT --- source/module_hsolver/test/hsolver_supplementary_mock.h | 6 ++++++ source/module_hsolver/test/test_hsolver_pw.cpp | 2 ++ source/module_hsolver/test/test_hsolver_sdft.cpp | 2 ++ 3 files changed, 10 insertions(+) diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index 95c69e555e..fd7d482d09 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -57,6 +57,12 @@ void ElecState::init_ks(Charge* chg_in, // pointer for class Charge return; } +Potential::~Potential(){} + +void Potential::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff){} + +void Potential::cal_fixed_v(double* vl_pseudo){} + } // namespace elecstate diff --git a/source/module_hsolver/test/test_hsolver_pw.cpp b/source/module_hsolver/test/test_hsolver_pw.cpp index cc545d8d14..5bbbc4bc0f 100644 --- a/source/module_hsolver/test/test_hsolver_pw.cpp +++ b/source/module_hsolver/test/test_hsolver_pw.cpp @@ -84,6 +84,7 @@ class TestHSolverPW : public ::testing::Test { TEST_F(TestHSolverPW, solve) { // initial memory and data elecstate_test.ekb.create(1, 2); + elecstate_test.pot = new elecstate::Potential; this->ekb_f.resize(2); psi_test_cf.resize(1, 2, 3); psi_test_cd.resize(1, 2, 3); @@ -216,6 +217,7 @@ TEST_F(TestHSolverPW, SolveLcaoInPW) { pwbk.nks = 1; // initial memory and data elecstate_test.ekb.create(1, 2); + elecstate_test.pot = new elecstate::Potential; // 1 kpt, 2 bands, 3 basis psi_test_cf.resize(1, 2, 3); psi_test_cd.resize(1, 2, 3); diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp index c48db2cee0..4e6aa18bd4 100644 --- a/source/module_hsolver/test/test_hsolver_sdft.cpp +++ b/source/module_hsolver/test/test_hsolver_sdft.cpp @@ -170,6 +170,7 @@ TEST_F(TestHSolverPW_SDFT, solve) { //initial memory and data elecstate_test.ekb.create(1,2); + elecstate_test.pot = new elecstate::Potential; elecstate_test.f_en.eband = 0.0; stowf.nbands_diag = 0; stowf.nbands_total = 0; @@ -213,6 +214,7 @@ TEST_F(TestHSolverPW_SDFT, solve_noband_skipcharge) { //initial memory and data elecstate_test.ekb.create(1,2); + elecstate_test.pot = new elecstate::Potential; elecstate_test.f_en.eband = 0.0; stowf.nbands_diag = 0; stowf.nbands_total = 0;