From 7a2e4a8dbaa21b794a1a581613b95ebc763a3de6 Mon Sep 17 00:00:00 2001 From: liiutao <74701833+A-006@users.noreply.github.com> Date: Wed, 18 Sep 2024 17:25:01 +0800 Subject: [PATCH 1/8] Refactor:Use PARAM instead of GlobalV::domag* (#5115) * add domag in systyem_parameter * change paramter noncolin in cpp file * change paramter noncolin in test file * change paramter lspinorb in cpp file * change paramter lspinorb in test file * change paramter npol in test file * change paramter npol in test file * change paramter domag_z in cpp file * change paramter domag_z in test file * change paramter domag in cpp file * change paramter domag in test file * fix bug * add unit test * delete parameter in variable * delete parameter out_pot in variable * delete nspin condition --------- Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> --- source/module_base/global_variable.cpp | 11 --- source/module_base/global_variable.h | 10 --- source/module_base/module_device/device.h | 2 +- .../module_basis/module_nao/beta_radials.cpp | 7 +- source/module_cell/atom_pseudo.cpp | 3 +- source/module_cell/read_atoms.cpp | 16 +--- source/module_cell/read_pp.cpp | 9 ++- source/module_cell/test/atom_pseudo_test.cpp | 2 +- source/module_cell/test/read_pp_test.cpp | 8 +- source/module_cell/test/unitcell_test.cpp | 4 +- .../module_cell/test/unitcell_test_para.cpp | 2 +- .../module_cell/test/unitcell_test_readpp.cpp | 8 +- source/module_cell/unitcell.cpp | 10 +-- source/module_elecstate/elecstate.cpp | 2 +- source/module_elecstate/elecstate_pw.cpp | 14 ++-- .../module_elecstate/kernels/elecstate_op.h | 5 +- .../module_elecstate/module_charge/charge.cpp | 14 ++-- .../module_charge/charge_mixing.cpp | 77 ++++++++++++------- .../potentials/potential_new.cpp | 6 +- .../test/charge_mixing_test.cpp | 20 ++--- .../test/elecstate_base_test.cpp | 4 +- .../test/elecstate_energy_test.cpp | 4 +- .../test/elecstate_print_test.cpp | 2 +- .../test/elecstate_pw_test.cpp | 8 +- .../test/potential_new_test.cpp | 30 +++++--- source/module_esolver/esolver.cpp | 8 +- source/module_esolver/esolver_fp.cpp | 4 +- source/module_esolver/esolver_ks.cpp | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 9 --- source/module_esolver/esolver_ks_pw.cpp | 2 +- source/module_esolver/lcao_before_scf.cpp | 4 +- source/module_esolver/pw_init_globalc.cpp | 2 +- .../module_xc/test/test_xc3.cpp | 5 +- .../module_xc/xc_functional_gradcorr.cpp | 11 +-- .../module_xc/xc_functional_vxc.cpp | 6 +- .../hamilt_lcaodft/LCAO_nl_mu.cpp | 2 +- .../hamilt_lcaodft/LCAO_set_st.cpp | 2 +- .../hamilt_lcaodft/fvnl_dbeta_k.cpp | 3 +- .../operator_lcao/op_exx_lcao.hpp | 4 +- .../hamilt_lcaodft/record_adj.cpp | 4 +- .../hamilt_lcaodft/spar_dh.cpp | 5 +- .../hamilt_lcaodft/spar_st.cpp | 5 +- .../hamilt_lcaodft/wavefunc_in_pw.cpp | 6 +- .../module_deepks/LCAO_deepks_pdm.cpp | 8 +- .../module_deepks/LCAO_deepks_psialpha.cpp | 10 +-- .../module_deepks/LCAO_deepks_torch.cpp | 4 +- .../module_deepks/cal_gdmx.cpp | 5 +- .../module_deepks/cal_gdmx_k.cpp | 5 +- .../module_deepks/deepks_fgamma.cpp | 4 +- .../module_deepks/deepks_fk.cpp | 4 +- .../module_deepks/orbital_precalc.cpp | 4 +- .../module_deepks/orbital_precalc_k.cpp | 4 +- .../module_deepks/test/nnr.cpp | 8 +- .../module_deepks/v_delta_precalc.cpp | 8 +- .../module_hamilt_lcao/module_dftu/dftu.cpp | 10 +-- .../module_dftu/dftu_folding.cpp | 12 +-- .../module_dftu/dftu_force.cpp | 4 +- .../module_dftu/dftu_io.cpp | 12 +-- .../module_dftu/dftu_occup.cpp | 13 ++-- .../module_dftu/dftu_tools.cpp | 9 ++- .../module_hamilt_lcao/module_gint/gint.cpp | 2 +- .../module_gint/gint_k_env.cpp | 5 +- .../module_gint/gint_k_pvpr.cpp | 7 +- .../module_gint/gint_k_sparse1.cpp | 13 ++-- .../module_gint/grid_technique.cpp | 6 +- .../hamilt_pwdft/VNL_in_pw.cpp | 26 +++---- .../module_hamilt_pw/hamilt_pwdft/elecond.cpp | 2 +- .../hamilt_pwdft/hamilt_pw.cpp | 2 +- .../hamilt_pwdft/structure_factor.cpp | 4 +- .../hamilt_pwdft/wavefunc.cpp | 22 +++--- .../hamilt_pwdft/wf_atomic.cpp | 14 ++-- .../hamilt_stodft/sto_forces.cpp | 4 +- .../hamilt_stodft/sto_hchi.cpp | 15 ++-- .../hamilt_stodft/sto_stress_pw.cpp | 21 +++-- source/module_hsolver/diago_elpa_native.cpp | 3 +- source/module_io/cal_r_overlap_R.cpp | 16 ++-- source/module_io/input_conv.cpp | 34 +------- source/module_io/read_set_globalv.cpp | 49 ++++++++++++ source/module_io/read_wfc_pw.cpp | 13 ++-- source/module_io/test/write_orb_info_test.cpp | 2 +- .../test_serial/io_system_variable_test.cpp | 16 ++++ source/module_io/to_wannier90_lcao.cpp | 6 +- source/module_io/to_wannier90_lcao_in_pw.cpp | 9 ++- source/module_io/unk_overlap_pw.cpp | 15 ++-- source/module_lr/potentials/pot_hxc_lrtd.cpp | 3 +- source/module_parameter/system_parameter.h | 7 +- source/module_psi/psi.cpp | 9 ++- source/module_psi/psi_initializer.cpp | 24 +++--- source/module_psi/psi_initializer_atomic.cpp | 2 +- source/module_psi/psi_initializer_nao.cpp | 2 +- .../test/psi_initializer_unit_test.cpp | 32 ++++---- 91 files changed, 460 insertions(+), 411 deletions(-) mode change 100755 => 100644 source/module_elecstate/module_charge/charge_mixing.cpp diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index a510e6fb87..21c636242e 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -73,22 +73,11 @@ std::ofstream ofs_info; // output math lib info std::ofstream ofs_device; // output device info -// added by zhengdy-soc -bool NONCOLIN = false; -bool LSPINORB = false; -bool DOMAG = false; -bool DOMAG_Z = false; -int NPOL = 1; - -std::vector rpa_orbitals; - //========================================================== // device flags added by denghui //========================================================== std::string device_flag = "unknown"; -int out_pot = 0; - double nelec = 0; diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 064198ca52..e582478ae8 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -35,12 +35,6 @@ extern bool use_uspp; extern std::string KS_SOLVER; // xiaohui add 2013-09-01 extern double SEARCH_RADIUS; // 11.1 // mohan add 2011-03-10 -// added by zhengdy-soc -extern bool NONCOLIN; // 0 : collinear ; 1 : non-collinear -extern bool LSPINORB; // 0 : no soc ; 1 : has soc -extern bool DOMAG; // 1 : calculate the magnetism with x, y, z component -extern bool DOMAG_Z; // 1 : constrain the magnetism to z axis -extern int NPOL; // 1 : no soc; 2 : has soc extern int PW_DIAG_NDIM; // 14 extern double PW_DIAG_THR; // 15 pw_diag_thr @@ -107,8 +101,6 @@ extern std::ofstream ofs_warning; extern std::ofstream ofs_info; extern std::ofstream ofs_device; -// rpa related -extern std::vector rpa_orbitals; // mixing parameters @@ -120,8 +112,6 @@ extern std::string device_flag; // precision flags added by denghui //========================================================== -extern int out_pot; - // "out_chg" elec step. /// @brief method to initialize wavefunction /// @author kirk0830, 20230920 diff --git a/source/module_base/module_device/device.h b/source/module_base/module_device/device.h index cb6e716e61..d0ed332555 100644 --- a/source/module_base/module_device/device.h +++ b/source/module_base/module_device/device.h @@ -44,7 +44,7 @@ int get_device_kpar(const int& kpar); /** * @brief Get the device flag object - * for module_io GlobalV::device_flag + * for module_io PARAM.globalv.device_flag */ std::string get_device_flag(const std::string& device, const std::string& ks_solver, diff --git a/source/module_basis/module_nao/beta_radials.cpp b/source/module_basis/module_nao/beta_radials.cpp index 63daa5cc09..fdef1822d0 100644 --- a/source/module_basis/module_nao/beta_radials.cpp +++ b/source/module_basis/module_nao/beta_radials.cpp @@ -1,5 +1,6 @@ #include "module_basis/module_nao/beta_radials.h" +#include "module_parameter/parameter.h" #include "module_base/global_variable.h" #include "module_base/parallel_common.h" #include "module_base/tool_quit.h" @@ -402,7 +403,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst //#endif // // // It is an error if lspinorb is set to true but the pseudopotential file does not contain spin-orbit information -// if (!has_so && GlobalV::LSPINORB) +// if (!has_so && PARAM.inp.lspinorb) // { // ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201", // "lspinorb is set to true but the pseudopotential file does not contain spin-orbit information"); @@ -558,7 +559,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst // double* rbeta_final = nullptr; // if (rank == 0) // { -// if (!GlobalV::LSPINORB && has_so) +// if (!PARAM.inp.lspinorb && has_so) // { // /* // * read the PP_DIJ block for averaging @@ -727,7 +728,7 @@ void BetaRadials::build(const Numerical_Nonlocal& nl, const int itype, std::ofst // delete[] ngrid_final; // delete[] rbeta_final; // -// if (rank == 0 && !GlobalV::LSPINORB && has_so) +// if (rank == 0 && !PARAM.inp.lspinorb && has_so) // { // delete[] l; // delete[] ngrid; diff --git a/source/module_cell/atom_pseudo.cpp b/source/module_cell/atom_pseudo.cpp index 7f15bce71e..0e4739fb77 100644 --- a/source/module_cell/atom_pseudo.cpp +++ b/source/module_cell/atom_pseudo.cpp @@ -1,5 +1,6 @@ #include "atom_pseudo.h" +#include "module_parameter/parameter.h" Atom_pseudo::Atom_pseudo() { } @@ -67,7 +68,7 @@ void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in, if (this->lmax > -1) { - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { int is = 0; for (int is1 = 0; is1 < 2; is1++) diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index bd09a0df54..2ebe9c29fb 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -139,18 +139,6 @@ int UnitCell::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running) } } - if (PARAM.globalv.rpa_setorb) - { - if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "ABFS_ORBITAL")) - { - std::cout << "RPA_EXX_LCAO read abfs_orb!!!" << std::endl; - GlobalV::rpa_orbitals.resize(ntype); - for (int i = 0; i < ntype; i++) - { - ifa >> GlobalV::rpa_orbitals[i]; - } - } - } #endif // __EXX #endif // __MPI #endif // __LCAO @@ -640,11 +628,9 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn if(GlobalV::NSPIN==4) { - if(GlobalV::NONCOLIN) + if(PARAM.inp.noncolin) { //if magnetization only along z-axis, default settings are DOMAG_Z=true and DOMAG=false - GlobalV::DOMAG_Z = false; - GlobalV::DOMAG = true; if(input_angle_mag) { atoms[it].m_loc_[ia].z = atoms[it].mag[ia] * diff --git a/source/module_cell/read_pp.cpp b/source/module_cell/read_pp.cpp index f314be2ec7..29c75afbe8 100644 --- a/source/module_cell/read_pp.cpp +++ b/source/module_cell/read_pp.cpp @@ -1,5 +1,6 @@ #include "read_pp.h" +#include "module_parameter/parameter.h" #include #include // Peize Lin fix bug about strcpy 2016-08-02 @@ -123,7 +124,7 @@ int Pseudopot_upf::average_p(const double& lambda, Atom_pseudo& pp) { int error = 0; double lambda_ = lambda; - if(!GlobalV::LSPINORB) { lambda_ = 0.0; } + if(!PARAM.inp.lspinorb) { lambda_ = 0.0; } if (pp.has_so && pp.tvanp) { error++; @@ -132,7 +133,7 @@ int Pseudopot_upf::average_p(const double& lambda, Atom_pseudo& pp) std::cout << "------------------------------------------------------" << std::endl; return error; } - if (!pp.has_so && GlobalV::LSPINORB) + if (!pp.has_so && PARAM.inp.lspinorb) { error++; std::cout << "warning_quit! no soc upf used for lspinorb calculation, error!" << std::endl; @@ -140,13 +141,13 @@ int Pseudopot_upf::average_p(const double& lambda, Atom_pseudo& pp) } // ModuleBase::WARNING_QUIT("average_p", "no soc upf used for lspinorb calculation, error!"); - if (!pp.has_so || (GlobalV::LSPINORB && std::abs(lambda_ - 1.0) < 1.0e-8)) + if (!pp.has_so || (PARAM.inp.lspinorb && std::abs(lambda_ - 1.0) < 1.0e-8)) { return error; } //if(std::abs(lambda_)<1.0e-8) - if(!GlobalV::LSPINORB) + if(!PARAM.inp.lspinorb) { int new_nbeta = 0; //calculate the new nbeta for(int nb=0; nb< pp.nbeta; nb++) diff --git a/source/module_cell/test/atom_pseudo_test.cpp b/source/module_cell/test/atom_pseudo_test.cpp index b93b2dcdf1..d1efd31295 100644 --- a/source/module_cell/test/atom_pseudo_test.cpp +++ b/source/module_cell/test/atom_pseudo_test.cpp @@ -58,7 +58,7 @@ TEST_F(AtomPseudoTest, SetDSo) atom_pseudo->set_d_so(d_so_in,nproj,nproj_soc,has_so); EXPECT_NEAR(atom_pseudo->d_so(0,0,0).real(),1e-8,1e-7); EXPECT_NEAR(atom_pseudo->d_so(0,0,0).imag(),1e-8,1e-7); - GlobalV::LSPINORB = true; + PARAM.input.lspinorb = true; atom_pseudo->set_d_so(d_so_in,nproj,nproj_soc,has_so); EXPECT_NEAR(atom_pseudo->d_so(0,0,0).real(),1e-8,1e-7); EXPECT_NEAR(atom_pseudo->d_so(0,0,0).imag(),1e-8,1e-7); diff --git a/source/module_cell/test/read_pp_test.cpp b/source/module_cell/test/read_pp_test.cpp index 86cf571fb2..bad80000f0 100644 --- a/source/module_cell/test/read_pp_test.cpp +++ b/source/module_cell/test/read_pp_test.cpp @@ -744,7 +744,7 @@ TEST_F(ReadPPTest, AverageSimpleReturns) int ierr; double lambda = 1.0; // first return - GlobalV::LSPINORB = 1; + PARAM.input.lspinorb = 1; upf->has_so = 0; ierr = read_pp->average_p(lambda, *upf); EXPECT_EQ(ierr,1); @@ -767,7 +767,7 @@ TEST_F(ReadPPTest, AverageErrReturns) ifs.open("./support/Te.pbe-rrkj.UPF"); read_pp->read_pseudo_upf(ifs, *upf); EXPECT_TRUE(upf->has_so); // has soc info - GlobalV::LSPINORB = 0; + PARAM.input.lspinorb = 0; ierr = read_pp->average_p(lambda, *upf); EXPECT_EQ(upf->nbeta,3); EXPECT_EQ(ierr,1); @@ -787,7 +787,7 @@ TEST_F(ReadPPTest, AverageLSPINORB0) int ierr; double lambda = 1.0; // LSPINORB = 0 - GlobalV::LSPINORB = 0; + PARAM.input.lspinorb = 0; ierr = read_pp->average_p(lambda, *upf); EXPECT_EQ(ierr,0); EXPECT_EQ(upf->nbeta,4); @@ -804,7 +804,7 @@ TEST_F(ReadPPTest, AverageLSPINORB1) int ierr; double lambda = 1.1; // LSPINORB = 0 - GlobalV::LSPINORB = 1; + PARAM.input.lspinorb = 1; ierr = read_pp->average_p(lambda, *upf); EXPECT_EQ(ierr,0); EXPECT_EQ(upf->nbeta,6); diff --git a/source/module_cell/test/unitcell_test.cpp b/source/module_cell/test/unitcell_test.cpp index 470590ac0e..c8d1509ce9 100644 --- a/source/module_cell/test/unitcell_test.cpp +++ b/source/module_cell/test/unitcell_test.cpp @@ -1320,7 +1320,7 @@ TEST_F(UcellTest, ReadAtomPositionsS4Noncolin) PARAM.input.basis_type = "lcao"; PARAM.sys.deepks_setorb = true; GlobalV::NSPIN = 4; - GlobalV::NONCOLIN = true; + PARAM.input.noncolin = true; EXPECT_NO_THROW(ucell->read_atom_species(ifa, ofs_running)); EXPECT_DOUBLE_EQ(ucell->latvec.e11, 4.27957); EXPECT_DOUBLE_EQ(ucell->latvec.e22, 4.27957); @@ -1352,7 +1352,7 @@ TEST_F(UcellTest, ReadAtomPositionsS4Colin) PARAM.input.basis_type = "lcao"; PARAM.sys.deepks_setorb = true; GlobalV::NSPIN = 4; - GlobalV::NONCOLIN = false; + PARAM.input.noncolin = false; EXPECT_NO_THROW(ucell->read_atom_species(ifa, ofs_running)); EXPECT_DOUBLE_EQ(ucell->latvec.e11, 4.27957); EXPECT_DOUBLE_EQ(ucell->latvec.e22, 4.27957); diff --git a/source/module_cell/test/unitcell_test_para.cpp b/source/module_cell/test/unitcell_test_para.cpp index 70a21ff5c8..76b76a3894 100644 --- a/source/module_cell/test/unitcell_test_para.cpp +++ b/source/module_cell/test/unitcell_test_para.cpp @@ -76,7 +76,7 @@ class UcellTest : public ::testing::Test PARAM.input.relax_new = utp.relax_new; PARAM.sys.global_out_dir = "./"; ucell = utp.SetUcellInfo(); - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; pp_dir = "./support/"; PARAM.input.pseudo_rcut = 15.0; PARAM.input.dft_functional = "default"; diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index 0713096d47..782ac08521 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/source/module_cell/test/unitcell_test_readpp.cpp @@ -110,7 +110,7 @@ class UcellTest : public ::testing::Test { PARAM.input.relax_new = utp.relax_new; PARAM.sys.global_out_dir = "./"; ucell = utp.SetUcellInfo(); - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; pp_dir = "./support/"; PARAM.input.pseudo_rcut = 15.0; PARAM.input.dft_functional = "default"; @@ -124,7 +124,7 @@ class UcellTest : public ::testing::Test { using UcellDeathTest = UcellTest; TEST_F(UcellDeathTest, ReadCellPPWarning1) { - GlobalV::LSPINORB = true; + PARAM.input.lspinorb = true; ucell->pseudo_fn[1] = "H_sr.upf"; testing::internal::CaptureStdout(); EXPECT_EXIT(ucell->read_cell_pseudopots(pp_dir, ofs), @@ -224,7 +224,7 @@ TEST_F(UcellTest, CalNatomwfc1) { } TEST_F(UcellTest, CalNatomwfc2) { - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; GlobalV::NSPIN = 4; ucell->read_cell_pseudopots(pp_dir, ofs); EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); @@ -238,7 +238,7 @@ TEST_F(UcellTest, CalNatomwfc2) { } TEST_F(UcellTest, CalNatomwfc3) { - GlobalV::LSPINORB = true; + PARAM.input.lspinorb = true; GlobalV::NSPIN = 4; ucell->read_cell_pseudopots(pp_dir, ofs); EXPECT_TRUE(ucell->atoms[0].ncpp.has_so); diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 068bb35c5e..bd36b66949 100644 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -579,10 +579,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { #ifdef __MPI Parallel_Common::bcast_bool(ok); Parallel_Common::bcast_bool(ok2); - if (GlobalV::NSPIN == 4) { - Parallel_Common::bcast_bool(GlobalV::DOMAG); - Parallel_Common::bcast_bool(GlobalV::DOMAG_Z); - } #endif if (!ok) { ModuleBase::WARNING_QUIT( @@ -985,15 +981,15 @@ void UnitCell::cal_nwfc(std::ofstream& log) { this->iwt2iw = new int[GlobalV::NLOCAL]; this->itia2iat.create(ntype, namax); - // this->itiaiw2iwt.create(ntype, namax, nwmax*GlobalV::NPOL); - this->set_iat2iwt(GlobalV::NPOL); + // this->itiaiw2iwt.create(ntype, namax, nwmax*PARAM.globalv.npol); + this->set_iat2iwt(PARAM.globalv.npol); int iat = 0; int iwt = 0; for (int it = 0; it < ntype; it++) { for (int ia = 0; ia < atoms[it].na; ia++) { this->itia2iat(it, ia) = iat; // this->iat2ia[iat] = ia; - for (int iw = 0; iw < atoms[it].nw * GlobalV::NPOL; iw++) { + for (int iw = 0; iw < atoms[it].nw * PARAM.globalv.npol; iw++) { // this->itiaiw2iwt(it, ia, iw) = iwt; this->iwt2iat[iwt] = iat; this->iwt2iw[iwt] = iw; diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index 8077abd43a..b2e0efc4d1 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -265,7 +265,7 @@ void ElecState::cal_nbands() // calculate number of bands (setup.f90) //======================================= double occupied_bands = static_cast(GlobalV::nelec / ModuleBase::DEGSPIN); - if (GlobalV::LSPINORB == 1) { + if (PARAM.inp.lspinorb == 1) { occupied_bands = static_cast(GlobalV::nelec); } diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index 155122c89f..78c699f1f3 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -47,7 +47,7 @@ ElecStatePW::~ElecStatePW() template void ElecStatePW::init_rho_data() { - if (GlobalV::device_flag == "gpu" || PARAM.inp.precision == "single") { + if (PARAM.globalv.device_flag == "gpu" || PARAM.inp.precision == "single") { this->rho = new Real*[this->charge->nspin]; resmem_var_op()(this->ctx, this->rho_data, this->charge->nspin * this->charge->nrxx); for (int ii = 0; ii < this->charge->nspin; ii++) { @@ -108,7 +108,7 @@ void ElecStatePW::psiToRho(const psi::Psi& psi) { this->add_usrho(psi); } - if (GlobalV::device_flag == "gpu" || PARAM.inp.precision == "single") { + if (PARAM.globalv.device_flag == "gpu" || PARAM.inp.precision == "single") { for (int ii = 0; ii < GlobalV::NSPIN; ii++) { castmem_var_d2h_op()(cpu_ctx, this->ctx, this->charge->rho[ii], this->rho[ii], this->charge->nrxx); if (get_xc_func_type() == 3) @@ -186,8 +186,8 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) { // replaced by denghui at 20221110 elecstate_pw_op()(this->ctx, - GlobalV::DOMAG, - GlobalV::DOMAG_Z, + PARAM.globalv.domag, + PARAM.globalv.domag_z, this->basis->nrxx, w1, this->rho, @@ -324,7 +324,7 @@ void ElecStatePW::add_usrho(const psi::Psi& psi) for (int ia = 0; ia < atom->na; ia++) { const int iat = ucell->itia2iat(it, ia); - if (GlobalV::NONCOLIN) + if (PARAM.inp.noncolin) { // noncolinear case } @@ -360,10 +360,10 @@ void ElecStatePW::add_usrho(const psi::Psi& psi) } // copy output from GEMM into desired format - if (GlobalV::NONCOLIN && !atom->ncpp.has_so) + if (PARAM.inp.noncolin && !atom->ncpp.has_so) { } - else if (GlobalV::NONCOLIN && atom->ncpp.has_so) + else if (PARAM.inp.noncolin && atom->ncpp.has_so) { } else diff --git a/source/module_elecstate/kernels/elecstate_op.h b/source/module_elecstate/kernels/elecstate_op.h index c1e1550f23..69c8cd33ea 100644 --- a/source/module_elecstate/kernels/elecstate_op.h +++ b/source/module_elecstate/kernels/elecstate_op.h @@ -1,5 +1,6 @@ // TODO: This is a temperary location for these functions. // And will be moved to a global module(module base) later. +#include "module_parameter/parameter.h" #ifndef MODULE_ELECSTATE_ELECSTATE_MULTI_DEVICE_H #define MODULE_ELECSTATE_ELECSTATE_MULTI_DEVICE_H #include @@ -32,8 +33,8 @@ struct elecstate_pw_op { /// /// Input Parameters /// @param ctx - which device this function runs on - /// @param DOMAG - GlobalV::DOMAG - /// @param DOMAG_Z - GlobalV::DOMAG_Z + /// @param DOMAG - PARAM.globalv.domag + /// @param DOMAG_Z - PARAM.globalv.domag_z /// @param nrxx - number of planewaves /// @param weight - input constant /// @param wfcr - input array, psi in real space diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index 66cf3733c4..78f653a762 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -539,7 +539,7 @@ void Charge::atomic_rho(const int spin_number_need, if (startmag_type == 1) { double sin_a1, sin_a2, cos_a1, cos_a2; - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { // will not be used now, will be deleted later ModuleBase::libm::sincos(atom->angle1[0], &sin_a1, &cos_a1); ModuleBase::libm::sincos(atom->angle2[0], &sin_a2, &cos_a2); @@ -551,7 +551,7 @@ void Charge::atomic_rho(const int spin_number_need, { const std::complex swap = strucFac(it, ig) * rho_lgl[this->rhopw->ig2igg[ig]]; rho_g3d(0, ig) += swap; - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { // will not be used now, will be deleted later rho_g3d(1, ig) += swap * (ucell.magnet.start_magnetization[it] / atom->ncpp.zv) * sin_a1 * cos_a2; @@ -560,7 +560,7 @@ void Charge::atomic_rho(const int spin_number_need, rho_g3d(3, ig) += swap * (ucell.magnet.start_magnetization[it] / atom->ncpp.zv) * cos_a1; } - else if (GlobalV::DOMAG_Z) + else if (PARAM.globalv.domag_z) { rho_g3d(1, ig) = 0.0; rho_g3d(2, ig) = 0.0; @@ -574,11 +574,11 @@ void Charge::atomic_rho(const int spin_number_need, for (int ia = 0; ia < atom->na; ia++) { double sin_a1, sin_a2, cos_a1, cos_a2; - if (GlobalV::DOMAG || GlobalV::DOMAG_Z) + if (PARAM.globalv.domag || PARAM.globalv.domag_z) { ModuleBase::libm::sincos(atom->angle1[ia], &sin_a1, &cos_a1); } - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { ModuleBase::libm::sincos(atom->angle2[ia], &sin_a2, &cos_a2); } @@ -596,12 +596,12 @@ void Charge::atomic_rho(const int spin_number_need, // calculate rho_total rho_g3d(0, ig) += swap; // calculate mag_z - if (GlobalV::DOMAG || GlobalV::DOMAG_Z) + if (PARAM.globalv.domag || PARAM.globalv.domag_z) { rho_g3d(3, ig) += swap * (atom->mag[ia] / atom->ncpp.zv) * cos_a1; } // calculate mag_x and mag_y - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { rho_g3d(1, ig) += swap * (atom->mag[ia] / atom->ncpp.zv) * sin_a1 * cos_a2; rho_g3d(2, ig) += swap * (atom->mag[ia] / atom->ncpp.zv) * sin_a1 * sin_a2; diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp old mode 100755 new mode 100644 index 38a0784693..cf6e6df4d5 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -161,7 +161,8 @@ void Charge_Mixing::init_mixing() // initailize nhat_mdata #ifdef USE_PAW - if(PARAM.inp.use_paw) this->mixing->init_mixing_data(this->nhat_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double)); + if(PARAM.inp.use_paw) { this->mixing->init_mixing_data(this->nhat_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double)); +} #endif ModuleBase::timer::tick("Charge_Mixing", "init_mixing"); @@ -239,7 +240,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) // The inner_product_real function (L1-norm) is different from that (L2-norm) in mixing. for (int is = 0; is < GlobalV::NSPIN; is++) { - if (is != 0 && is != 3 && GlobalV::DOMAG_Z) + if (is != 0 && is != 3 && PARAM.globalv.domag_z) { continue; } @@ -278,7 +279,7 @@ double Charge_Mixing::get_dkin(Charge* chr, const double nelec) // Get dkin from kin_r and kin_r_save for PW and LCAO both, which is different from drho. for (int is = 0; is < GlobalV::NSPIN; is++) { - if (is != 0 && is != 3 && GlobalV::DOMAG_Z) + if (is != 0 && is != 3 && PARAM.globalv.domag_z) { continue; } @@ -505,7 +506,8 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) { chr->rhog[0][ig] = rhog_magabs[ig]; // rhog double norm = std::sqrt(chr->rho[1][ig] * chr->rho[1][ig] + chr->rho[2][ig] * chr->rho[2][ig] + chr->rho[3][ig] * chr->rho[3][ig]); - if (std::abs(norm) < 1e-10) continue; + if (std::abs(norm) < 1e-10) { continue; +} double rescale_tmp = rho_magabs[npw + ig] / norm; chr->rho[1][ig] *= rescale_tmp; chr->rho[2][ig] *= rescale_tmp; @@ -770,7 +772,8 @@ void Charge_Mixing::mix_rho_real(Charge* chr) { chr->rho[0][ir] = rho_magabs[ir]; // rho double norm = std::sqrt(chr->rho[1][ir] * chr->rho[1][ir] + chr->rho[2][ir] * chr->rho[2][ir] + chr->rho[3][ir] * chr->rho[3][ir]); - if (norm < 1e-10) continue; + if (norm < 1e-10) { continue; +} double rescale_tmp = rho_magabs[nrxx + ir] / norm; chr->rho[1][ir] *= rescale_tmp; chr->rho[2][ir] *= rescale_tmp; @@ -1018,7 +1021,7 @@ void Charge_Mixing::mix_rho(Charge* chr) std::vector rho123(GlobalV::NSPIN * nrxx); for (int is = 0; is < GlobalV::NSPIN; ++is) { - if (is == 0 || is == 3 || !GlobalV::DOMAG_Z) + if (is == 0 || is == 3 || !PARAM.globalv.domag_z) { double* rho123_is = rho123.data() + is * nrxx; #ifdef _OPENMP @@ -1078,7 +1081,7 @@ void Charge_Mixing::mix_rho(Charge* chr) // rho_save is the charge before mixing for (int is = 0; is < GlobalV::NSPIN; ++is) { - if (is == 0 || is == 3 || !GlobalV::DOMAG_Z) + if (is == 0 || is == 3 || !PARAM.globalv.domag_z) { double* rho123_is = rho123.data() + is * nrxx; #ifdef _OPENMP @@ -1122,21 +1125,24 @@ void Charge_Mixing::mix_rho(Charge* chr) } #endif - if (new_e_iteration) + if (new_e_iteration) { new_e_iteration = false; +} ModuleBase::timer::tick("Charge", "mix_rho"); return; } void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) { - if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) + if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) { return; +} double fac, gg0, amin; // consider a resize for mixing_angle int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2; + if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) { resize_tmp = 2; +} // implement Kerker for density and magnetization separately for (int is = 0; is < GlobalV::NSPIN / resize_tmp; ++is) @@ -1181,11 +1187,13 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) void Charge_Mixing::Kerker_screen_real(double* drhor) { - if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) + if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) { return; +} // consider a resize for mixing_angle int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2; + if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) { resize_tmp = 2; +} // std::vector> drhog(this->rhopw->npw * GlobalV::NSPIN / resize_tmp); std::vector drhor_filter(this->rhopw->nrxx * GlobalV::NSPIN / resize_tmp); @@ -1208,7 +1216,8 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) assert(is == 1); // make sure break works #endif double is_mag = GlobalV::NSPIN - 1; - if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) is_mag = 1; + if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) { is_mag = 1; +} for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) { drhog[is * this->rhopw->npw + ig] = 0; @@ -1282,8 +1291,9 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) + if (this->rhopw->gg[ig] < 1e-8) { continue; +} sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1303,8 +1313,9 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) + if (this->rhopw->gg[ig] < 1e-8) { continue; +} sum += (conj(rhog1[0][ig] + rhog1[1][ig]) * (rhog2[0][ig] + rhog2[1][ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1343,9 +1354,9 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c } case 4: // non-collinear spin, added by zhengdy - if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) + if (!PARAM.globalv.domag && !PARAM.globalv.domag_z) { sum += part_of_noncolin(); - else + } else { // another part with magnetization #ifdef _OPENMP @@ -1353,8 +1364,9 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) + if (ig == this->rhopw->ig_gge0) { continue; +} sum += (conj(rhog1[0][ig]) * rhog2[0][ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1375,8 +1387,9 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == ig0) + if (ig == ig0) { continue; +} sum += fac3 * ((conj(rhog1[1][ig]) * rhog2[1][ig]).real() + (conj(rhog1[2][ig]) * rhog2[2][ig]).real() + (conj(rhog1[3][ig]) * rhog2[3][ig]).real()); @@ -1405,7 +1418,8 @@ double Charge_Mixing::inner_product_recip_simple(std::complex* rho1, std double rnorm = 0.0; // consider a resize for mixing_angle int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2; + if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) { resize_tmp = 2; +} #ifdef _OPENMP #pragma omp parallel for reduction(+ : rnorm) #endif @@ -1443,8 +1457,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) + if (this->rhopw->gg[ig] < 1e-8) { continue; +} sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1463,8 +1478,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { - if (this->rhopw->gg[ig] < 1e-8) + if (this->rhopw->gg[ig] < 1e-8) { continue; +} sum += (conj(rhog1[ig]) * (rhog2[ig])).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1500,7 +1516,7 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s } else if (GlobalV::NSPIN==4) { - if (!GlobalV::DOMAG && !GlobalV::DOMAG_Z) + if (!PARAM.globalv.domag && !PARAM.globalv.domag_z) { sum += part_of_rho(); } @@ -1512,8 +1528,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) + if (ig == this->rhopw->ig_gge0) { continue; +} sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1534,8 +1551,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == ig0) + if (ig == ig0) { continue; +} sum += fac3 * ((conj(rhog1[ig + npw]) * rhog2[ig + npw]).real() + (conj(rhog1[ig + 2*npw]) * rhog2[ig + 2*npw]).real() + (conj(rhog1[ig + 3*npw]) * rhog2[ig + 3*npw]).real()); @@ -1549,8 +1567,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == this->rhopw->ig_gge0) + if (ig == this->rhopw->ig_gge0) { continue; +} sum += (conj(rhog1[ig]) * rhog2[ig]).real() / this->rhopw->gg[ig]; } sum *= fac; @@ -1570,8 +1589,9 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s #endif for (int ig = 0; ig < this->rhopw->npw; ig++) { - if (ig == ig0) + if (ig == ig0) { continue; +} sum += fac3 * ((conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real()); } @@ -1593,7 +1613,8 @@ double Charge_Mixing::inner_product_real(double* rho1, double* rho2) double rnorm = 0.0; // consider a resize for mixing_angle int resize_tmp = 1; - if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) resize_tmp = 2; + if (GlobalV::NSPIN == 4 && this->mixing_angle > 0) { resize_tmp = 2; +} #ifdef _OPENMP #pragma omp parallel for reduction(+ : rnorm) diff --git a/source/module_elecstate/potentials/potential_new.cpp b/source/module_elecstate/potentials/potential_new.cpp index da81b9ecbf..6c6270b24e 100644 --- a/source/module_elecstate/potentials/potential_new.cpp +++ b/source/module_elecstate/potentials/potential_new.cpp @@ -46,7 +46,7 @@ Potential::~Potential() } this->components.clear(); } - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { delmem_sd_op()(gpu_ctx, s_veff_smooth); delmem_sd_op()(gpu_ctx, s_vofk_smooth); @@ -129,7 +129,7 @@ void Potential::allocate() this->vofk_smooth.create(GlobalV::NSPIN, nrxx_smooth); ModuleBase::Memory::record("Pot::vofk_smooth", sizeof(double) * GlobalV::NSPIN * nrxx_smooth); } - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { resmem_sd_op()(gpu_ctx, s_veff_smooth, GlobalV::NSPIN * nrxx_smooth); resmem_sd_op()(gpu_ctx, s_vofk_smooth, GlobalV::NSPIN * nrxx_smooth); @@ -177,7 +177,7 @@ void Potential::update_from_charge(const Charge*const chg, const UnitCell*const } #endif - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { castmem_d2s_h2d_op()(gpu_ctx, cpu_ctx, diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index c0c3798369..210cf922c8 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -430,13 +430,13 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest) drhog2[i] = std::complex(1.0, 1.0); } - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); EXPECT_NEAR(inner, 28260.091995611871, 1e-8); PARAM.sys.gamma_only_pw= true; - GlobalV::DOMAG = true; - GlobalV::DOMAG_Z = true; + PARAM.sys.domag = true; + PARAM.sys.domag_z = true; inner = CMtest.inner_product_recip_hartree(drhog1.data(), drhog2.data()); EXPECT_NEAR(inner, 110668.61166927818, 1e-8); @@ -526,13 +526,13 @@ TEST_F(ChargeMixingTest, InnerDotRecipRhoTest) drhog2[i] = std::complex(1.0, 1.0); } - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); EXPECT_NEAR(inner, 28260.091995611871, 1e-8); PARAM.sys.gamma_only_pw= true; - GlobalV::DOMAG = true; - GlobalV::DOMAG_Z = true; + PARAM.sys.domag = true; + PARAM.sys.domag_z = true; inner = CMtest.inner_product_recip_rho(drhog1.data(), drhog2.data()); EXPECT_NEAR(inner, 110668.61166927818, 1e-8); } @@ -750,7 +750,7 @@ TEST_F(ChargeMixingTest, MixRhoTest) PARAM.sys.double_grid = false; charge.set_rhopw(&pw_basis); const int nspin = GlobalV::NSPIN = 1; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag_z = false; FUNC_TYPE = 3; PARAM.input.mixing_beta = 0.7; PARAM.input.mixing_ndim = 1; @@ -884,7 +884,7 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) PARAM.sys.double_grid = true; charge.set_rhopw(&pw_dbasis); const int nspin = GlobalV::NSPIN = 1; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag_z = false; FUNC_TYPE = 3; PARAM.input.mixing_beta = 0.7; PARAM.input.mixing_ndim = 1; diff --git a/source/module_elecstate/test/elecstate_base_test.cpp b/source/module_elecstate/test/elecstate_base_test.cpp index dd2dff8a0d..984b0a1777 100644 --- a/source/module_elecstate/test/elecstate_base_test.cpp +++ b/source/module_elecstate/test/elecstate_base_test.cpp @@ -134,7 +134,7 @@ class MockElecState : public ElecState GlobalV::NBANDS = 6; GlobalV::NLOCAL = 6; PARAM.input.esolver_type = "ksdft"; - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; PARAM.input.basis_type = "pw"; GlobalV::KPAR = 1; GlobalV::NPROC_IN_POOL = 1; @@ -198,7 +198,7 @@ TEST_F(ElecStateTest, CalNbandsFractionElec) TEST_F(ElecStateTest, CalNbandsSOC) { - GlobalV::LSPINORB = true; + PARAM.input.lspinorb = true; GlobalV::NBANDS = 0; elecstate->cal_nbands(); EXPECT_EQ(GlobalV::NBANDS, 20); diff --git a/source/module_elecstate/test/elecstate_energy_test.cpp b/source/module_elecstate/test/elecstate_energy_test.cpp index 5c7fd676f3..200220aba7 100644 --- a/source/module_elecstate/test/elecstate_energy_test.cpp +++ b/source/module_elecstate/test/elecstate_energy_test.cpp @@ -83,11 +83,11 @@ class MockElecState : public ElecState GlobalV::NBANDS = 6; GlobalV::NLOCAL = 6; PARAM.input.esolver_type = "ksdft"; - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; PARAM.input.basis_type = "pw"; GlobalV::KPAR = 1; GlobalV::NPROC_IN_POOL = 1; - PARAM.input.sc_mag_switch = 1; + PARAM.input.sc_mag_switch = true; } }; const double* ElecState::getRho(int spin) const diff --git a/source/module_elecstate/test/elecstate_print_test.cpp b/source/module_elecstate/test/elecstate_print_test.cpp index ec1c006998..44317d4376 100644 --- a/source/module_elecstate/test/elecstate_print_test.cpp +++ b/source/module_elecstate/test/elecstate_print_test.cpp @@ -387,7 +387,7 @@ TEST_F(ElecStatePrintTest, PrintEtotColorS4) GlobalV::TWO_EFERMI = true; PARAM.input.out_bandgap = true; GlobalV::NSPIN = 4; - GlobalV::NONCOLIN = true; + PARAM.input.noncolin = true; GlobalV::MY_RANK = 0; elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); } diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp index e31908f1b5..e059518fc8 100644 --- a/source/module_elecstate/test/elecstate_pw_test.cpp +++ b/source/module_elecstate/test/elecstate_pw_test.cpp @@ -164,10 +164,10 @@ void Charge::check_rho() void Set_GlobalV_Default() { - GlobalV::device_flag = "cpu"; + PARAM.sys.device_flag = "cpu"; PARAM.input.precision = "double"; - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; // Base class dependent GlobalV::NSPIN = 1; GlobalV::nelec = 10.0; @@ -176,7 +176,7 @@ void Set_GlobalV_Default() GlobalV::NBANDS = 6; GlobalV::NLOCAL = 6; PARAM.input.esolver_type = "ksdft"; - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; PARAM.input.basis_type = "pw"; GlobalV::KPAR = 1; GlobalV::NPROC_IN_POOL = 1; diff --git a/source/module_elecstate/test/potential_new_test.cpp b/source/module_elecstate/test/potential_new_test.cpp index c53ff6a58e..8468f79213 100644 --- a/source/module_elecstate/test/potential_new_test.cpp +++ b/source/module_elecstate/test/potential_new_test.cpp @@ -55,7 +55,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) void Set_GlobalV_Default() { GlobalV::NSPIN = 1; - GlobalV::device_flag = "cpu"; + PARAM.sys.device_flag = "cpu"; PARAM.input.precision = "double"; } } // namespace elecstate @@ -117,22 +117,30 @@ class PotentialNewTest : public ::testing::Test } virtual void TearDown() { - if (rhopw != nullptr) + if (rhopw != nullptr) { delete rhopw; - if (rhodpw != nullptr) +} + if (rhodpw != nullptr) { delete rhodpw; - if (ucell != nullptr) +} + if (ucell != nullptr) { delete ucell; - if (vloc != nullptr) +} + if (vloc != nullptr) { delete vloc; - if (structure_factors != nullptr) +} + if (structure_factors != nullptr) { delete structure_factors; - if (etxc != nullptr) +} + if (etxc != nullptr) { delete etxc; - if (vtxc != nullptr) +} + if (vtxc != nullptr) { delete vtxc; - if (pot != nullptr) +} + if (pot != nullptr) { delete pot; +} } }; @@ -185,7 +193,7 @@ TEST_F(PotentialNewTest, ConstructorGPUDouble) { // this is just a trivial call to the GPU code rhopw->nrxx = 100; - GlobalV::device_flag = "gpu"; + PARAM.sys.device_flag = "gpu"; pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc); EXPECT_TRUE(pot->fixed_mode); EXPECT_TRUE(pot->dynamic_mode); @@ -198,7 +206,7 @@ TEST_F(PotentialNewTest, ConstructorGPUSingle) { // this is just a trivial call to the GPU code rhopw->nrxx = 100; - GlobalV::device_flag = "gpu"; + PARAM.sys.device_flag = "gpu"; PARAM.input.precision = "single"; pot = new elecstate::Potential(rhopw, rhopw, ucell, vloc, structure_factors, etxc, vtxc); EXPECT_TRUE(pot->fixed_mode); diff --git a/source/module_esolver/esolver.cpp b/source/module_esolver/esolver.cpp index f5bfc64f0c..67f1b5d95e 100644 --- a/source/module_esolver/esolver.cpp +++ b/source/module_esolver/esolver.cpp @@ -100,7 +100,7 @@ std::string determine_type() GlobalV::ofs_running << " The esolver type has been set to : " << esolver_type << std::endl; - auto device_info = GlobalV::device_flag; + auto device_info = PARAM.globalv.device_flag; for (char &c : device_info) { @@ -112,11 +112,11 @@ std::string determine_type() if (GlobalV::MY_RANK == 0) { std::cout << " RUNNING WITH DEVICE : " << device_info << " / " - << base_device::information::get_device_info(GlobalV::device_flag) << std::endl; + << base_device::information::get_device_info(PARAM.globalv.device_flag) << std::endl; } GlobalV::ofs_running << "\n RUNNING WITH DEVICE : " << device_info << " / " - << base_device::information::get_device_info(GlobalV::device_flag) << std::endl; + << base_device::information::get_device_info(PARAM.globalv.device_flag) << std::endl; return esolver_type; } @@ -132,7 +132,7 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) if (esolver_type == "ksdft_pw") { #if ((defined __CUDA) || (defined __ROCM)) - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index f9104b8d76..331192ed85 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -15,11 +15,11 @@ namespace ModuleESolver ESolver_FP::ESolver_FP() { // pw_rho = new ModuleBase::PW_Basis(); - pw_rho = new ModulePW::PW_Basis_Big(GlobalV::device_flag, PARAM.inp.precision); + pw_rho = new ModulePW::PW_Basis_Big(PARAM.globalv.device_flag, PARAM.inp.precision); if ( PARAM.globalv.double_grid) { - pw_rhod = new ModulePW::PW_Basis_Big(GlobalV::device_flag, PARAM.inp.precision); + pw_rhod = new ModulePW::PW_Basis_Big(PARAM.globalv.device_flag, PARAM.inp.precision); } else { diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index cf3457cc62..e71661b00c 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -55,7 +55,7 @@ ESolver_KS::ESolver_KS() // pw_rho = new ModuleBase::PW_Basis(); // temporary, it will be removed - pw_wfc = new ModulePW::PW_Basis_K_Big(GlobalV::device_flag, PARAM.inp.precision); + pw_wfc = new ModulePW::PW_Basis_K_Big(PARAM.globalv.device_flag, PARAM.inp.precision); ModulePW::PW_Basis_K_Big* tmp = static_cast(pw_wfc); // should not use INPUT here, mohan 2024-05-12 diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 0b4f96a18a..56c564dd97 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -880,14 +880,6 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) istep); } - // 3) print potential - if (this->conv_elec || iter == PARAM.inp.scf_nmax) - { - if (GlobalV::out_pot < 0) // mohan add 2011-10-10 - { - GlobalV::out_pot = -2; - } - } if (!this->conv_elec) { @@ -1173,7 +1165,6 @@ void ESolver_KS_LCAO::after_scf(const int istep) { // ModuleRPA::DFT_RPA_interface // rpa_interface(GlobalC::exx_info.info_global); - // rpa_interface.rpa_exx_lcao().info.files_abfs = GlobalV::rpa_orbitals; RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), MPI_COMM_WORLD, diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 556f281632..a5dbacac83 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -213,7 +213,7 @@ void ESolver_KS_PW::before_scf(const int istep) //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho - if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) + if (GlobalV::NSPIN == 4 && PARAM.globalv.domag) { GlobalC::ucell.cal_ux(); } diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index e071ac33b5..685880699e 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -152,7 +152,7 @@ void ESolver_KS_LCAO::beforesolver(const int istep) PARAM.inp.sc_mag_switch, GlobalC::ucell, PARAM.inp.sc_file, - GlobalV::NPOL, + PARAM.globalv.npol, &(this->pv), GlobalV::NSPIN, this->kv, @@ -165,7 +165,7 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho //========================================================= - if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) + if (GlobalV::NSPIN == 4 && PARAM.globalv.domag) { GlobalC::ucell.cal_ux(); } diff --git a/source/module_esolver/pw_init_globalc.cpp b/source/module_esolver/pw_init_globalc.cpp index c74cbebc13..9fd7b7110e 100644 --- a/source/module_esolver/pw_init_globalc.cpp +++ b/source/module_esolver/pw_init_globalc.cpp @@ -90,7 +90,7 @@ void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, &(this->sf)); this->kspw_psi - = GlobalV::device_flag == "gpu" || PARAM.inp.precision == "single" + = PARAM.globalv.device_flag == "gpu" || PARAM.inp.precision == "single" ? new psi::Psi(this->psi[0]) : reinterpret_cast*>(this->psi); diff --git a/source/module_hamilt_general/module_xc/test/test_xc3.cpp b/source/module_hamilt_general/module_xc/test/test_xc3.cpp index e9d3bd59b2..70cfba4d9b 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc3.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc3.cpp @@ -1,5 +1,8 @@ #include "gtest/gtest.h" #include "xctest.h" +#define private public +#include "module_parameter/parameter.h" +#undef private #include "../xc_functional.h" #include "../exx_info.h" #include "xc3_mock.h" @@ -90,7 +93,7 @@ class XCTest_GRADCORR : public XCTest XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true); GlobalV::NSPIN = 4; - GlobalV::DOMAG = true; + PARAM.sys.domag = true; XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false); } }; diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp index 4d94f68ec4..68cf4efb3b 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp @@ -1,5 +1,6 @@ // This file contains subroutines realted to gradient calculations // it contains 5 subroutines: +#include "module_parameter/parameter.h" // 1. gradcorr, which calculates gradient correction // 2. grad_wfc, which calculates gradient of wavefunction // it is used in stress_func_mgga.cpp @@ -34,7 +35,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, int nspin0 = GlobalV::NSPIN; if(GlobalV::NSPIN==4) { nspin0 =1; } - if(GlobalV::NSPIN==4&&(GlobalV::DOMAG||GlobalV::DOMAG_Z)) { nspin0 = 2; + if(GlobalV::NSPIN==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) { nspin0 = 2; } assert(nspin0>0); @@ -124,7 +125,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); } - if(GlobalV::NSPIN == 4&&(GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if(GlobalV::NSPIN == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) { rhotmp2 = new double[rhopw->nrxx]; rhogsum2 = new std::complex[rhopw->npw]; @@ -389,7 +390,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, else { double zeta = ( rhotmp1[ir] - rhotmp2[ir] ) / rh; - if(GlobalV::NSPIN==4&&(GlobalV::DOMAG||GlobalV::DOMAG_Z)) { zeta = fabs(zeta) * neg[ir]; + if(GlobalV::NSPIN==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) { zeta = fabs(zeta) * neg[ir]; } const double grh2 = (gdr1[ir]+gdr2[ir]).norm2(); XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); @@ -546,7 +547,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, vtxc += vtxcgc; etxc += etxcgc; - if(GlobalV::NSPIN == 4 && (GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if(GlobalV::NSPIN == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) { #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 1024) @@ -591,7 +592,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, if(!is_stress) { delete[] h2; } } - if(GlobalV::NSPIN == 4 && (GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if(GlobalV::NSPIN == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) { delete[] neg; if(!is_stress) diff --git a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp index 3f377dd2df..67cafccd23 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp @@ -37,7 +37,7 @@ std::tuple XC_Functional::v_xc( double vanishing_charge = 1.0e-10; - if (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) + if (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) { // spin-unpolarized case #ifdef _OPENMP @@ -188,7 +188,7 @@ std::tuple XC_Functional::v_xc_libxc( // Peiz ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); const int nspin = - (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) + (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) ? 1 : 2; double etxc = 0.0; @@ -450,7 +450,7 @@ std::tuple XC_Functional::v_xc_libxc( // Peiz for( int ir=0; ir < beta | phi2 > diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp index c4188bce32..cc28a075b0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_set_st.cpp @@ -324,7 +324,7 @@ void build_ST_new(ForceStressArrays& fsr, ModuleBase::timer::tick("LCAO_domain", "build_ST_new"); const int nspin = GlobalV::NSPIN; - const int npol = GlobalV::NPOL; + const int npol = PARAM.globalv.npol; const bool cal_stress = GlobalV::CAL_STRESS; const bool gamma_only_local = PARAM.globalv.gamma_only_local; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp index dc694e4eb5..ce7c407446 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_k.cpp @@ -1,5 +1,6 @@ #include "FORCE.h" #include "module_base/memory.h" +#include "module_parameter/parameter.h" #include "module_base/parallel_reduce.h" #include "module_base/timer.h" #include "module_base/tool_threading.h" @@ -41,7 +42,7 @@ void Force_LCAO>::cal_fvnl_dbeta(const elecstate::DensityMa const int nspin = GlobalV::NSPIN; const int nspin_DMR = (nspin == 2) ? 2 : 1; - const int npol = GlobalV::NPOL; + const int npol = PARAM.globalv.npol; // Data structure for storing , for a detailed description // check out the same data structure in build_Nonlocal_mu_new diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp index 58a8c13a69..1b3ad2bc34 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp @@ -212,7 +212,7 @@ void OperatorEXX>::contributeHR() GlobalC::exx_info.info_global.hybrid_alpha, *this->Hexxd, *this->hR->get_paraV(), - GlobalV::NPOL, + PARAM.globalv.npol, *this->hR, this->use_cell_nearest ? &this->cell_nearest : nullptr); } else { @@ -221,7 +221,7 @@ void OperatorEXX>::contributeHR() GlobalC::exx_info.info_global.hybrid_alpha, *this->Hexxc, *this->hR->get_paraV(), - GlobalV::NPOL, + PARAM.globalv.npol, *this->hR, this->use_cell_nearest ? &this->cell_nearest : nullptr); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp index 03729d09a8..43569e2e24 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp @@ -127,7 +127,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto ++na_each[iat]; if (!gamma_only) { - for(int ii=0; iinw * GlobalV::NPOL; ++ii) + for(int ii=0; iinw * PARAM.globalv.npol; ++ii) { // the index of orbitals in this processor const int iw1_all = start1 + ii; @@ -135,7 +135,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto if(mu<0) {continue; } - for(int jj=0; jj @@ -178,7 +179,7 @@ void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, Abfs::Vector3_Order dR(grid.getBox(ad).x, grid.getBox(ad).y, grid.getBox(ad).z); - for (int ii = 0; ii < atom1->nw * GlobalV::NPOL; ii++) + for (int ii = 0; ii < atom1->nw * PARAM.globalv.npol; ii++) { const int iw1_all = start + ii; const int mu = pv.global2local_row(iw1_all); @@ -188,7 +189,7 @@ void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, continue; } - for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++) + for (int jj = 0; jj < atom2->nw * PARAM.globalv.npol; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp index 686689499b..8070eb2d53 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp @@ -1,5 +1,6 @@ #include "spar_st.h" +#include "module_parameter/parameter.h" #include "force_stress_arrays.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" // only for INPUT @@ -158,7 +159,7 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, grid.getBox(ad).y, grid.getBox(ad).z); - for (int ii = 0; ii < atom1->nw * GlobalV::NPOL; ii++) { + for (int ii = 0; ii < atom1->nw * PARAM.globalv.npol; ii++) { const int iw1_all = start + ii; const int mu = pv.global2local_row(iw1_all); @@ -166,7 +167,7 @@ void sparse_format::cal_STN_R_for_T(const UnitCell& ucell, continue; } - for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++) { + for (int jj = 0; jj < atom2->nw * PARAM.globalv.npol; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp index e471cb4ade..e5776de524 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp @@ -301,7 +301,7 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, if(GlobalC::ucell.atoms[it].ncpp.has_so) { const double j = std::abs(double(L+is_N) - 0.5); - if (!(GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if (!(PARAM.globalv.domag||PARAM.globalv.domag_z)) {//atomic_wfc_so for(int m=0; m<2*L+1; m++) { @@ -388,7 +388,7 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, iwall++; } iwall += 2*L +1; - } // end else INPUT.starting_spin_angle || !GlobalV::DOMAG + } // end else INPUT.starting_spin_angle || !PARAM.globalv.domag } // end if GlobalC::ucell.atoms[it].has_so else {//atomic_wfc_nc @@ -432,7 +432,7 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, iwall += 2*L+1; } // end else GlobalC::ucell.atoms[it].has_so } // end for is_N - } // end if GlobalV::NONCOLIN + } // end if PARAM.inp.noncolin else {//LSDA and nomagnet case for(int m=0; m<2*L+1; m++) diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp index 1983df875e..2925dfc69a 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp @@ -167,7 +167,7 @@ void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); const double dist1 = (tau1-tau0).norm() * ucell.lat0; if (dist1 >= Rcut_Alpha + Rcut_AO1) @@ -200,7 +200,7 @@ void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; const double Rcut_AO2 = orb.Phi[T2].getRcut(); const double dist2 = (tau2-tau0).norm() * ucell.lat0; @@ -417,7 +417,7 @@ void LCAO_Deepks::cal_projected_DM_k(const elecstate::DensityMatrix tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); const double dist1 = (tau1-tau0).norm() * ucell.lat0; if (dist1 >= Rcut_Alpha + Rcut_AO1) @@ -453,7 +453,7 @@ void LCAO_Deepks::cal_projected_DM_k(const elecstate::DensityMatrix tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; ModuleBase::Vector3 dR2(GridD.getBox(ad2).x, GridD.getBox(ad2).y, GridD.getBox(ad2).z); const double Rcut_AO2 = orb.Phi[T2].getRcut(); diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_psialpha.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_psialpha.cpp index 83bf92dd59..393456c4f1 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_psialpha.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_psialpha.cpp @@ -85,9 +85,9 @@ void LCAO_Deepks::build_psialpha(const bool& calc_deri, all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end()); // middle loop : all atomic basis on the adjacent atom ad - for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += GlobalV::NPOL) + for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += PARAM.globalv.npol) { - const int iw1 = all_indexes[iw1l] / GlobalV::NPOL; + const int iw1 = all_indexes[iw1l] / PARAM.globalv.npol; std::vector> nlm; // 2D, dim 0 contains the overlap // dim 1-3 contains the gradient of overlap @@ -109,7 +109,7 @@ void LCAO_Deepks::build_psialpha(const bool& calc_deri, else { nlm_cur.insert({all_indexes[iw1l], nlm}); - if (GlobalV::NPOL == 2) + if (PARAM.globalv.npol == 2) nlm_cur.insert({all_indexes[iw1l + 1], nlm}); } } // end iw @@ -189,7 +189,7 @@ void LCAO_Deepks::check_psialpha(const bool& calc_deri, const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw * GlobalV::NPOL; + const int nw1_tot = atom1->nw * PARAM.globalv.npol; const double dist1 = (tau1 - tau0).norm() * ucell.lat0; @@ -229,7 +229,7 @@ void LCAO_Deepks::check_psialpha(const bool& calc_deri, const int iw2_local = pv->global2local_col(iw1_all); if (iw1_local < 0 && iw2_local < 0) continue; - const int iw1_0 = iw1 / GlobalV::NPOL; + const int iw1_0 = iw1 / PARAM.globalv.npol; std::vector> nlm; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_torch.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_torch.cpp index 43b8426d76..cadc3a79f1 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_torch.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_torch.cpp @@ -187,7 +187,7 @@ void LCAO_Deepks::prepare_psialpha(const int nlocal, const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double dist1 = (tau1-tau0).norm() * ucell.lat0; @@ -204,7 +204,7 @@ void LCAO_Deepks::prepare_psialpha(const int nlocal, const int iw2_local = pv->global2local_col(iw1_all); if(iw1_local < 0 || iw2_local < 0) {continue; } - const int iw1_0 = iw1/GlobalV::NPOL; + const int iw1_0 = iw1/PARAM.globalv.npol; std::vector nlm = this->nlm_save[iat][ad][iw1][0]; int ib=0; diff --git a/source/module_hamilt_lcao/module_deepks/cal_gdmx.cpp b/source/module_hamilt_lcao/module_deepks/cal_gdmx.cpp index fe462bb0d2..3ecefaeada 100644 --- a/source/module_hamilt_lcao/module_deepks/cal_gdmx.cpp +++ b/source/module_hamilt_lcao/module_deepks/cal_gdmx.cpp @@ -1,5 +1,6 @@ #ifdef __DEEPKS +#include "module_parameter/parameter.h" #include "LCAO_deepks.h" #include "module_base/vector3.h" #include "module_base/timer.h" @@ -66,7 +67,7 @@ void LCAO_Deepks::cal_gdmx(const std::vector& dm, const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); for (int ad2=0; ad2 < GridD.getAdjacentNum()+1 ; ad2++) @@ -77,7 +78,7 @@ void LCAO_Deepks::cal_gdmx(const std::vector& dm, const int ibt2 = ucell.itia2iat(T2,I2); const ModuleBase::Vector3 tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; const double Rcut_AO2 = orb.Phi[T2].getRcut(); const double dist1 = (tau1-tau0).norm() * ucell.lat0; diff --git a/source/module_hamilt_lcao/module_deepks/cal_gdmx_k.cpp b/source/module_hamilt_lcao/module_deepks/cal_gdmx_k.cpp index 2baf9bdb92..eed8597534 100644 --- a/source/module_hamilt_lcao/module_deepks/cal_gdmx_k.cpp +++ b/source/module_hamilt_lcao/module_deepks/cal_gdmx_k.cpp @@ -1,5 +1,6 @@ #ifdef __DEEPKS +#include "module_parameter/parameter.h" #include "LCAO_deepks.h" #include "module_base/vector3.h" #include "module_base/timer.h" @@ -61,7 +62,7 @@ void LCAO_Deepks::cal_gdmx_k(const std::vector> const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); ModuleBase::Vector3 dR1(GridD.getBox(ad1).x, GridD.getBox(ad1).y, GridD.getBox(ad1).z); @@ -74,7 +75,7 @@ void LCAO_Deepks::cal_gdmx_k(const std::vector> const int ibt2 = ucell.itia2iat(T2,I2); const ModuleBase::Vector3 tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; ModuleBase::Vector3 dR2(GridD.getBox(ad2).x, GridD.getBox(ad2).y, GridD.getBox(ad2).z); const double Rcut_AO2 = orb.Phi[T2].getRcut(); diff --git a/source/module_hamilt_lcao/module_deepks/deepks_fgamma.cpp b/source/module_hamilt_lcao/module_deepks/deepks_fgamma.cpp index 8b190e2571..cebd0f7c21 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_fgamma.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_fgamma.cpp @@ -56,7 +56,7 @@ void DeePKS_domain::cal_f_delta_gamma( const int ibt1 = ucell.itia2iat(T1,I1); const ModuleBase::Vector3 tau1 = gd.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); for (int ad2=0; ad2 < gd.getAdjacentNum()+1 ; ad2++) @@ -66,7 +66,7 @@ void DeePKS_domain::cal_f_delta_gamma( const int ibt2 = ucell.itia2iat(T2,I2); const ModuleBase::Vector3 tau2 = gd.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; const double Rcut_AO2 = orb.Phi[T2].getRcut(); const double dist1 = (tau1-tau0).norm() * ucell.lat0; diff --git a/source/module_hamilt_lcao/module_deepks/deepks_fk.cpp b/source/module_hamilt_lcao/module_deepks/deepks_fk.cpp index 8d6413a331..26491736bf 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_fk.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_fk.cpp @@ -56,7 +56,7 @@ void DeePKS_domain::cal_f_delta_k( const int start1 = ucell.itiaiw2iwt(T1, I1, 0); const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); ModuleBase::Vector3 dR1(GridD.getBox(ad1).x, GridD.getBox(ad1).y, GridD.getBox(ad1).z); @@ -69,7 +69,7 @@ void DeePKS_domain::cal_f_delta_k( const int start2 = ucell.itiaiw2iwt(T2, I2, 0); const ModuleBase::Vector3 tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; ModuleBase::Vector3 dR2(GridD.getBox(ad2).x, GridD.getBox(ad2).y, GridD.getBox(ad2).z); const double Rcut_AO2 = orb.Phi[T2].getRcut(); diff --git a/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp b/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp index ec3b5d84ea..7a4b403605 100644 --- a/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp +++ b/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp @@ -76,7 +76,7 @@ void LCAO_Deepks::cal_orbital_precalc( const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw * GlobalV::NPOL; + const int nw1_tot = atom1->nw * PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); const double dist1 = (tau1 - tau0).norm() * ucell.lat0; @@ -116,7 +116,7 @@ void LCAO_Deepks::cal_orbital_precalc( const ModuleBase::Vector3 tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw * GlobalV::NPOL; + const int nw2_tot = atom2->nw * PARAM.globalv.npol; const double Rcut_AO2 = orb.Phi[T2].getRcut(); const double dist2 = (tau2 - tau0).norm() * ucell.lat0; diff --git a/source/module_hamilt_lcao/module_deepks/orbital_precalc_k.cpp b/source/module_hamilt_lcao/module_deepks/orbital_precalc_k.cpp index b632c02494..956ed15ba4 100644 --- a/source/module_hamilt_lcao/module_deepks/orbital_precalc_k.cpp +++ b/source/module_hamilt_lcao/module_deepks/orbital_precalc_k.cpp @@ -76,7 +76,7 @@ void LCAO_Deepks::cal_orbital_precalc_k( } const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw * GlobalV::NPOL; + const int nw1_tot = atom1->nw * PARAM.globalv.npol; ModuleBase::Vector3 dR1(GridD.getBox(ad1).x, GridD.getBox(ad1).y, @@ -126,7 +126,7 @@ void LCAO_Deepks::cal_orbital_precalc_k( const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw * GlobalV::NPOL; + const int nw2_tot = atom2->nw * PARAM.globalv.npol; ModuleBase::Vector3 dR2(GridD.getBox(ad2).x, GridD.getBox(ad2).y, diff --git a/source/module_hamilt_lcao/module_deepks/test/nnr.cpp b/source/module_hamilt_lcao/module_deepks/test/nnr.cpp index 35b3eb2e2b..b1d882da72 100644 --- a/source/module_hamilt_lcao/module_deepks/test/nnr.cpp +++ b/source/module_hamilt_lcao/module_deepks/test/nnr.cpp @@ -1,5 +1,7 @@ #include "LCAO_deepks_test.h" - +#define private public +#include "module_parameter/parameter.h" +#undef private void test_deepks::cal_nnr(void) { ModuleBase::TITLE("test_deepks","cal_nnr"); @@ -198,14 +200,14 @@ void test_deepks::folding_nnr(const Test_Deepks::K_Vectors &kv) // calculate how many matrix elements are in // this processor. //-------------------------------------------------- - for(int ii=0; iinw*GlobalV::NPOL; ii++) + for(int ii=0; iinw*PARAM.sys.npol; ii++) { // the index of orbitals in this processor const int iw1_all = start + ii; const int mu = ParaO.global2local_row(iw1_all); if(mu<0)continue; - for(int jj=0; jjnw*GlobalV::NPOL; jj++) + for(int jj=0; jjnw*PARAM.sys.npol; jj++) { int iw2_all = start2 + jj; const int nu = ParaO.global2local_col(iw2_all); diff --git a/source/module_hamilt_lcao/module_deepks/v_delta_precalc.cpp b/source/module_hamilt_lcao/module_deepks/v_delta_precalc.cpp index 8595ea55b5..0174e4137c 100644 --- a/source/module_hamilt_lcao/module_deepks/v_delta_precalc.cpp +++ b/source/module_hamilt_lcao/module_deepks/v_delta_precalc.cpp @@ -49,7 +49,7 @@ void LCAO_Deepks::cal_v_delta_precalc(const int nlocal, const int start1 = ucell.itiaiw2iwt(T1, I1, 0); const ModuleBase::Vector3 tau1 = GridD.getAdjacentTau(ad1); const Atom* atom1 = &ucell.atoms[T1]; - const int nw1_tot = atom1->nw*GlobalV::NPOL; + const int nw1_tot = atom1->nw*PARAM.globalv.npol; const double Rcut_AO1 = orb.Phi[T1].getRcut(); const double dist1 = (tau1-tau0).norm() * ucell.lat0; @@ -65,7 +65,7 @@ void LCAO_Deepks::cal_v_delta_precalc(const int nlocal, const int start2 = ucell.itiaiw2iwt(T2, I2, 0); const ModuleBase::Vector3 tau2 = GridD.getAdjacentTau(ad2); const Atom* atom2 = &ucell.atoms[T2]; - const int nw2_tot = atom2->nw*GlobalV::NPOL; + const int nw2_tot = atom2->nw*PARAM.globalv.npol; const double Rcut_AO2 = orb.Phi[T2].getRcut(); const double dist2 = (tau2-tau0).norm() * ucell.lat0; @@ -81,14 +81,14 @@ void LCAO_Deepks::cal_v_delta_precalc(const int nlocal, const int iw1_local = pv->global2local_row(iw1_all); if(iw1_local < 0) {continue; } - const int iw1_0 = iw1/GlobalV::NPOL; + const int iw1_0 = iw1/PARAM.globalv.npol; for (int iw2=0; iw2global2local_col(iw2_all); if(iw2_local < 0) {continue; } - const int iw2_0 = iw2/GlobalV::NPOL; + const int iw2_0 = iw2/PARAM.globalv.npol; std::vector nlm1 = this->nlm_save[iat][ad1][iw1][0]; std::vector nlm2 = this->nlm_save[iat][ad2][iw2][0]; diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp index abe68ad417..f463110d1d 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu.cpp @@ -56,7 +56,7 @@ void DFTU::init(UnitCell& cell, // unitcell class // needs reconstructions in future // global parameters, need to be removed in future - const int npol = GlobalV::NPOL; // number of polarization directions + const int npol = PARAM.globalv.npol; // number of polarization directions const int nlocal = GlobalV::NLOCAL; // number of total local orbitals const int nspin = GlobalV::NSPIN; // number of spins @@ -288,14 +288,14 @@ void DFTU::cal_energy_correction(const int istep) for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < GlobalV::NPOL; ipol0++) + for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) { const int m0_all = m0 + (2 * l + 1) * ipol0; nm_trace += this->locale[iat][l][n][0](m0_all, m0_all); for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { int m1_all = m1 + (2 * l + 1) * ipol1; @@ -319,12 +319,12 @@ void DFTU::cal_energy_correction(const int istep) // calculate the double counting term included in eband for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { const int m1_all = m1 + ipol1 * (2 * l + 1); for (int m2 = 0; m2 < 2 * l + 1; m2++) { - for (int ipol2 = 0; ipol2 < GlobalV::NPOL; ipol2++) + for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) { const int m2_all = m2 + ipol2 * (2 * l + 1); diff --git a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp index 9de5ff795c..ae6e99573b 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp @@ -92,9 +92,9 @@ void DFTU::fold_dSR_gamma( if (adj) { - for (int jj = 0; jj < atom1->nw * GlobalV::NPOL; ++jj) + for (int jj = 0; jj < atom1->nw * PARAM.globalv.npol; ++jj) { - const int jj0 = jj / GlobalV::NPOL; + const int jj0 = jj / PARAM.globalv.npol; const int iw1_all = start1 + jj0; const int mu = pv.global2local_row(iw1_all); if (mu < 0) @@ -102,9 +102,9 @@ void DFTU::fold_dSR_gamma( continue; } - for (int kk = 0; kk < atom2->nw * GlobalV::NPOL; ++kk) + for (int kk = 0; kk < atom2->nw * PARAM.globalv.npol; ++kk) { - const int kk0 = kk / GlobalV::NPOL; + const int kk0 = kk / PARAM.globalv.npol; const int iw2_all = start2 + kk0; const int nu = pv.global2local_col(iw2_all); if (nu < 0) @@ -224,7 +224,7 @@ void DFTU::folding_matrix_k( // calculate how many matrix elements are in // this processor. //-------------------------------------------------- - for (int ii = 0; ii < atom1->nw * GlobalV::NPOL; ii++) + for (int ii = 0; ii < atom1->nw * PARAM.globalv.npol; ii++) { // the index of orbitals in this processor const int iw1_all = start1 + ii; @@ -232,7 +232,7 @@ void DFTU::folding_matrix_k( if (mu < 0) { continue; } - for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++) + for (int jj = 0; jj < atom2->nw * PARAM.globalv.npol; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); diff --git a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp index 2dd344c687..0e19dcc3b5 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp @@ -353,7 +353,7 @@ void DFTU::cal_force_k(ForceStressArrays& fsr, for (int m = 0; m < 2 * l + 1; m++) { - for (int ipol = 0; ipol < GlobalV::NPOL; ipol++) + for (int ipol = 0; ipol < PARAM.globalv.npol; ipol++) { const int iwt = this->iatlnmipol2iwt[iat][l][n][m][ipol]; const int mu = pv.global2local_row(iwt); @@ -566,7 +566,7 @@ void DFTU::cal_force_gamma(const double* rho_VU, // Calculate the local occupation number matrix for (int m = 0; m < 2 * l + 1; m++) { - for (int ipol = 0; ipol < GlobalV::NPOL; ipol++) + for (int ipol = 0; ipol < PARAM.globalv.npol; ipol++) { const int iwt = this->iatlnmipol2iwt[iat][l][n][m][ipol]; const int mu = pv.global2local_row(iwt); diff --git a/source/module_hamilt_lcao/module_dftu/dftu_io.cpp b/source/module_hamilt_lcao/module_dftu/dftu_io.cpp index 0b92ecc1e8..490e5fcbd6 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_io.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_io.cpp @@ -194,13 +194,13 @@ void DFTU::write_occup_m(std::ofstream &ofs, bool diag) else { for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < GlobalV::NPOL; ipol0++) + for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) { const int m0_all = m0 + (2 * l + 1) * ipol0; for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { int m1_all = m1 + (2 * l + 1) * ipol1; ofs << std::setw(12) << std::setprecision(8) << std::fixed @@ -335,13 +335,13 @@ void DFTU::read_occup_m(const std::string &fn) double value = 0.0; for (int m0 = 0; m0 < 2 * L + 1; m0++) { - for (int ipol0 = 0; ipol0 < GlobalV::NPOL; ipol0++) + for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) { const int m0_all = m0 + (2 * L + 1) * ipol0; for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { int m1_all = m1 + (2 * L + 1) * ipol1; ifdftu >> value; @@ -431,13 +431,13 @@ void DFTU::local_occup_bcast() { for (int m0 = 0; m0 < 2 * L + 1; m0++) { - for (int ipol0 = 0; ipol0 < GlobalV::NPOL; ipol0++) + for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) { const int m0_all = m0 + (2 * L + 1) * ipol0; for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { int m1_all = m1 + (2 * L + 1) * ipol1; #ifdef __MPI diff --git a/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp b/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp index ae710d93dd..656300e9a7 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp @@ -1,5 +1,6 @@ #include "dftu.h" #include "module_base/timer.h" +#include "module_parameter/parameter.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" @@ -223,7 +224,7 @@ void DFTU::cal_occup_m_k(const int iter, // Calculate the local occupation number matrix for (int m0 = 0; m0 < 2 * l + 1; m0++) { - for (int ipol0 = 0; ipol0 < GlobalV::NPOL; ipol0++) + for (int ipol0 = 0; ipol0 < PARAM.globalv.npol; ipol0++) { const int iwt0 = this->iatlnmipol2iwt[iat][l][n][m0][ipol0]; const int mu = this->paraV->global2local_row(iwt0); @@ -231,7 +232,7 @@ void DFTU::cal_occup_m_k(const int iter, for (int m1 = 0; m1 < 2 * l + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { const int iwt1 = this->iatlnmipol2iwt[iat][l][n][m1][ipol1]; const int nu = this->paraV->global2local_col(iwt1); @@ -291,7 +292,7 @@ void DFTU::cal_occup_m_k(const int iter, ModuleBase::matrix temp(locale[iat][l][n][0]); MPI_Allreduce(&temp(0, 0), &locale[iat][l][n][0](0, 0), - (2 * l + 1) * GlobalV::NPOL * (2 * l + 1) * GlobalV::NPOL, + (2 * l + 1) * PARAM.globalv.npol * (2 * l + 1) * PARAM.globalv.npol, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); @@ -419,7 +420,7 @@ void DFTU::cal_occup_m_gamma(const int iter, const std::vectoriatlnmipol2iwt[iat][l][n][m0][ipol0]; const int mu = this->paraV->global2local_row(iwt0); @@ -427,7 +428,7 @@ void DFTU::cal_occup_m_gamma(const int iter, const std::vectoriatlnmipol2iwt[iat][l][n][m1][ipol1]; const int nu = this->paraV->global2local_col(iwt1); @@ -461,7 +462,7 @@ void DFTU::cal_occup_m_gamma(const int iter, const std::vectorparaV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) { @@ -41,7 +42,7 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com for (int m2 = 0; m2 < 2 * L + 1; m2++) { - for (int ipol2 = 0; ipol2 < GlobalV::NPOL; ipol2++) + for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); @@ -93,7 +94,7 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) for (int m1 = 0; m1 < 2 * L + 1; m1++) { - for (int ipol1 = 0; ipol1 < GlobalV::NPOL; ipol1++) + for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++) { const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) { @@ -101,7 +102,7 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) } for (int m2 = 0; m2 < 2 * L + 1; m2++) { - for (int ipol2 = 0; ipol2 < GlobalV::NPOL; ipol2++) + for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); diff --git a/source/module_hamilt_lcao/module_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/gint.cpp index 3c4e4b59cd..ba2e814a23 100644 --- a/source/module_hamilt_lcao/module_gint/gint.cpp +++ b/source/module_hamilt_lcao/module_gint/gint.cpp @@ -43,7 +43,7 @@ void Gint::cal_gint(Gint_inout* inout) { } if (this->gridt->max_atom > 0) { #ifdef __CUDA - if (GlobalV::device_flag == "gpu" + if (PARAM.globalv.device_flag == "gpu" && (inout->job == Gint_Tools::job_type::vlocal || inout->job == Gint_Tools::job_type::rho || inout->job == Gint_Tools::job_type::force)) { diff --git a/source/module_hamilt_lcao/module_gint/gint_k_env.cpp b/source/module_hamilt_lcao/module_gint/gint_k_env.cpp index 0557131cf3..023bc71c39 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_env.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_env.cpp @@ -1,5 +1,6 @@ #include "gint_k.h" #include "grid_technique.h" +#include "module_parameter/parameter.h" #include "module_base/timer.h" #include "module_base/ylm.h" #include "module_basis/module_ao/ORB_read.h" @@ -113,8 +114,8 @@ void Gint_k::cal_env_k(int ik, { for (int is = 0; is < 2; ++is) { - iw1_lo = this->gridt->trace_lo[start1] / GlobalV::NPOL - + this->gridt->lgd / GlobalV::NPOL * is; + iw1_lo = this->gridt->trace_lo[start1] / PARAM.globalv.npol + + this->gridt->lgd / PARAM.globalv.npol * is; for (int iw = 0; iw < atom1->nw; ++iw, ++iw1_lo) { tmp += std::complex(psi1[iw], 0.0) * psi_k[iw1_lo] * kphase; diff --git a/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp b/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp index 12ec052356..d1aa6624c7 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp @@ -1,5 +1,6 @@ #include "gint_k.h" #include "grid_technique.h" +#include "module_parameter/parameter.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/libm/libm.h" @@ -70,7 +71,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer* hR, const UnitCell* ucell } this->hRGint->set_zero(); - const int npol = GlobalV::NPOL; + const int npol = PARAM.globalv.npol; const UnitCell& ucell = *ucell_in; for (int iat = 0; iat < ucell.nat; ++iat) { @@ -194,7 +195,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer>* hR, const U } this->hRGintCd->set_zero(); - const int npol = GlobalV::NPOL; + const int npol = PARAM.globalv.npol; const UnitCell& ucell = *ucell_in; for (int iat = 0; iat < ucell.nat; ++iat) @@ -280,7 +281,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer>* hR, const U } tmp_pointer += 2 * atom2->nw; } - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { tmp_pointer = tmp_matrix->get_pointer(); for (int iw = 0; iw < atom1->nw; iw++) diff --git a/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp b/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp index 49678592e0..0314f83c4f 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp @@ -1,5 +1,6 @@ #include "gint_k.h" #include "grid_technique.h" +#include "module_parameter/parameter.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/memory.h" @@ -378,13 +379,13 @@ void Gint_k::cal_dvlocal_R_sparseMatrix(const int& current_spin, int ixxx = DM_start + this->gridt->find_R2st[iat][nad2]; - for (int iw = 0; iw < atom1->nw * GlobalV::NPOL; iw++) + for (int iw = 0; iw < atom1->nw * PARAM.globalv.npol; iw++) { - for (int iw2 = 0; iw2 < atom2->nw * GlobalV::NPOL; iw2++) + for (int iw2 = 0; iw2 < atom2->nw * PARAM.globalv.npol; iw2++) { const int nw = atom2->nw; - const int mug0 = iw / GlobalV::NPOL; - const int nug0 = iw2 / GlobalV::NPOL; + const int mug0 = iw / PARAM.globalv.npol; + const int nug0 = iw2 / PARAM.globalv.npol; const int iw_nowg = ixxx + mug0 * nw + nug0; if (GlobalV::NSPIN == 4) @@ -454,7 +455,7 @@ void Gint_k::cal_dvlocal_R_sparseMatrix(const int& current_spin, else if (iw % 2 == 0 && iw2 % 2 == 1) { // spin = 1; - if (!GlobalV::DOMAG) + if (!PARAM.globalv.domag) { // do nothing } @@ -489,7 +490,7 @@ void Gint_k::cal_dvlocal_R_sparseMatrix(const int& current_spin, else if (iw % 2 == 1 && iw2 % 2 == 0) { // spin = 2; - if (!GlobalV::DOMAG) + if (!PARAM.globalv.domag) { // do nothing } diff --git a/source/module_hamilt_lcao/module_gint/grid_technique.cpp b/source/module_hamilt_lcao/module_gint/grid_technique.cpp index 8c544434ee..d1610639d2 100644 --- a/source/module_hamilt_lcao/module_gint/grid_technique.cpp +++ b/source/module_hamilt_lcao/module_gint/grid_technique.cpp @@ -13,7 +13,7 @@ Grid_Technique::Grid_Technique() { allocate_find_R2 = false; #if ((defined __CUDA) /* || (defined __ROCM) */) - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { is_malloced = false; } #endif @@ -22,7 +22,7 @@ Grid_Technique::Grid_Technique() { Grid_Technique::~Grid_Technique() { #if ((defined __CUDA) /* || (defined __ROCM) */) - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { free_gpu_gint_variables(this->nat); } #endif @@ -118,7 +118,7 @@ void Grid_Technique::set_pbc_grid( this->cal_trace_lo(ucell); #if ((defined __CUDA) /* || (defined __ROCM) */) - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { this->init_gpu_gint_variables(ucell, num_stream); } #endif diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp index c91f23ccf7..6c8ab20973 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp @@ -32,7 +32,7 @@ void pseudopot_cell_vnl::release_memory() if (this->nhm <= 0 || memory_released) { return; } - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { @@ -153,7 +153,7 @@ void pseudopot_cell_vnl::init(const int ntype, this->deeq_nc.create(GlobalV::NSPIN, GlobalC::ucell.nat, this->nhm, this->nhm); this->qq_nt.create(ntype, this->nhm, this->nhm); this->qq_so.create(ntype, 4, this->nhm, this->nhm); - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { @@ -272,7 +272,7 @@ void pseudopot_cell_vnl::init(const int ntype, ModuleBase::Memory::record("VNL::tab_at", ntype * nchix_nc * GlobalV::NQX * sizeof(double)); } } - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { @@ -466,7 +466,7 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex { _gk[ig] = this->wfcpw->getgpluskcar(ik, ig); } - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { resmem_int_op()(ctx, atom_nh, GlobalC::ucell.ntype); resmem_int_op()(ctx, atom_nb, GlobalC::ucell.ntype); @@ -575,7 +575,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ // In the spin-orbit case we need the unitary matrix u which rotates the // real spherical harmonics and yields the complex ones. soc.fcoef.create(cell.ntype, this->nhm, this->nhm); - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { soc.rot_ylm(this->lmaxkb); } @@ -702,7 +702,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ { const int ir = static_cast(indv(it, ip)); const int is = static_cast(indv(it, ip2)); - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { this->dvan_so(0, it, ip, ip2) = cell.atoms[it].ncpp.dion(ir, is); this->dvan_so(3, it, ip, ip2) = cell.atoms[it].ncpp.dion(ir, is); @@ -784,7 +784,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ for (int jh = ih; jh < upf->nh; jh++) { this->radial_fft_q(1, ih, jh, it, &qnorm, ylmk0, &qgm); - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { this->qq_so(it, 0, ih, jh) = cell.omega * qgm.real(); this->qq_so(it, 0, jh, ih) = this->qq_so(it, 0, ih, jh); @@ -863,7 +863,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ delete[] aux; delete[] jl; } - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { @@ -1402,7 +1402,7 @@ void pseudopot_cell_vnl::cal_effective_D(const ModuleBase::matrix& veff, { for (int jh = ih; jh < nht; jh++) { - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { this->deeq_nc(is, iat, ih, jh) = this->dvan_so(is, it, ih, jh); this->deeq_nc(is, iat, jh, ih) = this->dvan_so(is, it, jh, ih); @@ -1454,7 +1454,7 @@ void pseudopot_cell_vnl::cal_effective_D(const ModuleBase::matrix& veff, for (int iat = 0; iat < cell.nat; iat++) { int it = cell.iat2it[iat]; - if (GlobalV::NONCOLIN) + if (PARAM.inp.noncolin) { if (cell.atoms[it].ncpp.has_so) { @@ -1481,7 +1481,7 @@ void pseudopot_cell_vnl::cal_effective_D(const ModuleBase::matrix& veff, } } } - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { @@ -1670,7 +1670,7 @@ void pseudopot_cell_vnl::newd_so(const int& iat, UnitCell& cell) { for (int lh = 0; lh < upf->nh; lh++) { - if (GlobalV::DOMAG) + if (PARAM.globalv.domag) { deeq_nc(ijs, iat, ih, jh) += deeq(0, iat, kh, lh) @@ -1713,7 +1713,7 @@ void pseudopot_cell_vnl::newd_nc(const int& iat, UnitCell& cell) { for (int jh = 0; jh < upf->nh; jh++) { - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { deeq_nc(0, iat, ih, jh) = dvan_so(0, it, ih, jh) + deeq(0, iat, ih, jh) + deeq(3, iat, ih, jh); deeq_nc(3, iat, ih, jh) = dvan_so(3, it, ih, jh) + deeq(0, iat, ih, jh) - deeq(3, iat, ih, jh); diff --git a/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp b/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp index 5b8d47cea7..4307018260 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp @@ -114,7 +114,7 @@ void EleCond::jjresponse_ks(const int ik, const int nt, const double dt, const d std::vector> pij(nbands * nbands); std::vector pij2(reducenb2, 0); // px|right> - velop.act(this->p_psi, nbands * GlobalV::NPOL, levc, prevc.data()); + velop.act(this->p_psi, nbands * PARAM.globalv.npol, levc, prevc.data()); for (int id = 0; id < ndim; ++id) { diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp index 7c68a7af6a..725591c990 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp @@ -276,7 +276,7 @@ void HamiltPW::sPsi(const T* psi_in, // psi setmem_complex_op()(this->ctx, ps, 0, this->ppcell->nkb * nbands); // spsi = psi + sum qq |beta> - if (GlobalV::NONCOLIN) + if (PARAM.inp.noncolin) { // spsi_nc std::cout << " noncolinear in uspp is not implemented yet " << std::endl; diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp index 02dca3e50f..575cb53508 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp @@ -21,7 +21,7 @@ Structure_Factor::Structure_Factor() Structure_Factor::~Structure_Factor() { - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { delmem_cd_op()(gpu_ctx, this->c_eigts1); delmem_cd_op()(gpu_ctx, this->c_eigts2); @@ -145,7 +145,7 @@ void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::P inat++; } } - if (GlobalV::device_flag == "gpu") { + if (PARAM.globalv.device_flag == "gpu") { if (PARAM.inp.precision == "single") { resmem_cd_op()(gpu_ctx, this->c_eigts1, Ucell->nat * (2 * rho_basis->nx + 1)); resmem_cd_op()(gpu_ctx, this->c_eigts2, Ucell->nat * (2 * rho_basis->ny + 1)); diff --git a/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp b/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp index 6decf17403..96084eda1e 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp @@ -46,7 +46,7 @@ psi::Psi>* wavefunc::allocate(const int nkstot, const int n int prefactor = 1; if (GlobalV::NSPIN == 4) { - prefactor = GlobalV::NPOL; // added by zhengdy-soc + prefactor = PARAM.globalv.npol; // added by zhengdy-soc } const int nks2 = nks; @@ -55,15 +55,15 @@ psi::Psi>* wavefunc::allocate(const int nkstot, const int n if (PARAM.inp.calculation == "nscf" && this->mem_saver == 1) { // initial psi rather than evc - psi_out = new psi::Psi>(1, GlobalV::NBANDS, npwx * GlobalV::NPOL, ngk); + psi_out = new psi::Psi>(1, GlobalV::NBANDS, npwx * PARAM.globalv.npol, ngk); if (PARAM.inp.basis_type == "lcao_in_pw") { - wanf2[0].create(GlobalV::NLOCAL, npwx * GlobalV::NPOL); - const size_t memory_cost = GlobalV::NLOCAL * (GlobalV::NPOL * npwx) * sizeof(std::complex); + wanf2[0].create(GlobalV::NLOCAL, npwx * PARAM.globalv.npol); + const size_t memory_cost = GlobalV::NLOCAL * (PARAM.globalv.npol * npwx) * sizeof(std::complex); std::cout << " Memory for wanf2 (MB): " << double(memory_cost) / 1024.0 / 1024.0 << std::endl; ModuleBase::Memory::record("WF::wanf2", memory_cost); } - const size_t memory_cost = GlobalV::NBANDS * (GlobalV::NPOL * npwx) * sizeof(std::complex); + const size_t memory_cost = GlobalV::NBANDS * (PARAM.globalv.npol * npwx) * sizeof(std::complex); std::cout << " MEMORY FOR PSI (MB) : " << double(memory_cost) / 1024.0 / 1024.0 << std::endl; ModuleBase::Memory::record("Psi_PW", memory_cost); } @@ -79,10 +79,10 @@ psi::Psi>* wavefunc::allocate(const int nkstot, const int n for (int ik = 0; ik < nks2; ik++) { - this->wanf2[ik].create(GlobalV::NLOCAL, npwx * GlobalV::NPOL); + this->wanf2[ik].create(GlobalV::NLOCAL, npwx * PARAM.globalv.npol); } - const size_t memory_cost = nks2 * GlobalV::NLOCAL * (npwx * GlobalV::NPOL) * sizeof(std::complex); + const size_t memory_cost = nks2 * GlobalV::NLOCAL * (npwx * PARAM.globalv.npol) * sizeof(std::complex); std::cout << " Memory for wanf2 (MB): " << double(memory_cost) / 1024.0 / 1024.0 << std::endl; ModuleBase::Memory::record("WF::wanf2", memory_cost); } @@ -90,8 +90,8 @@ psi::Psi>* wavefunc::allocate(const int nkstot, const int n else { // initial psi rather than evc - psi_out = new psi::Psi>(nks2, GlobalV::NBANDS, npwx * GlobalV::NPOL, ngk); - const size_t memory_cost = nks2 * GlobalV::NBANDS * (GlobalV::NPOL * npwx) * sizeof(std::complex); + psi_out = new psi::Psi>(nks2, GlobalV::NBANDS, npwx * PARAM.globalv.npol, ngk); + const size_t memory_cost = nks2 * GlobalV::NBANDS * (PARAM.globalv.npol * npwx) * sizeof(std::complex); std::cout << " MEMORY FOR PSI (MB) : " << double(memory_cost) / 1024.0 / 1024.0 << std::endl; ModuleBase::Memory::record("Psi_PW", memory_cost); } @@ -117,7 +117,7 @@ void wavefunc::wfcinit(psi::Psi>* psi_in, ModulePW::PW_Basi this->irindex = new int[wfc_basis->fftnxy]; wfc_basis->getfftixy2is(this->irindex); #if defined(__CUDA) || defined(__ROCM) - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { wfc_basis->get_ig2ixyz_k(); } @@ -731,7 +731,7 @@ void wavefunc::init_after_vc(const int nks) assert(GlobalV::NBANDS > 0); const int nks2 = nks; - const int nbasis = this->npwx * GlobalV::NPOL; + const int nbasis = this->npwx * PARAM.globalv.npol; if ((PARAM.inp.basis_type == "lcao" || PARAM.inp.basis_type == "lcao_in_pw") || winput::out_spillage == 2) { diff --git a/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp b/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp index 21405127fd..35ca3b564d 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp @@ -321,7 +321,7 @@ void WF_atomic::atomic_wfc(const int ik, Soc soc; soc.rot_ylm(l+1); const double j = GlobalC::ucell.atoms[it].ncpp.jchi[iw]; - if ( !(GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if ( !(PARAM.globalv.domag||PARAM.globalv.domag_z)) {//atomic_wfc_so double fact[2]; for(int m=-l-1;m* psi, FPTYPE *tmparg = new FPTYPE[nstnz]; for (int iw = iw_start ;iw < iw_end;iw++) { - std::complex* ppsi = &(psi[iw * this->npwx * GlobalV::NPOL]); + std::complex* ppsi = &(psi[iw * this->npwx * PARAM.globalv.npol]); int startig = 0; - for(int ipol = 0 ; ipol < GlobalV::NPOL ; ++ipol) + for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) { for(int ir=0; ir < nxy; ir++) @@ -662,7 +662,7 @@ void WF_atomic::random_t(std::complex* psi, #endif // __MPI for (int iw = iw_start ;iw < iw_end;iw++) { - std::complex* ppsi = &(psi[iw * this->npwx * GlobalV::NPOL]); + std::complex* ppsi = &(psi[iw * this->npwx * PARAM.globalv.npol]); for (int ig = 0;ig < ng;ig++) { const FPTYPE rr = std::rand()/FPTYPE(RAND_MAX); //qianrui add RAND_MAX @@ -670,7 +670,7 @@ void WF_atomic::random_t(std::complex* psi, const FPTYPE gk2 = wfc_basis->getgk2(ik,ig); ppsi[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); } - if(GlobalV::NPOL==2) {for (int ig = this->npwx;ig < this->npwx + ng;ig++) + if(PARAM.globalv.npol==2) {for (int ig = this->npwx;ig < this->npwx + ng;ig++) { const FPTYPE rr = std::rand()/FPTYPE(RAND_MAX); const FPTYPE arg= ModuleBase::TWO_PI * std::rand()/FPTYPE(RAND_MAX); @@ -710,7 +710,7 @@ void WF_atomic::atomicrandom(ModuleBase::ComplexMatrix& psi, for (int iw = iw_start ;iw < iw_end;iw++) { int startig = 0; - for(int ipol = 0 ; ipol < GlobalV::NPOL ; ++ipol) + for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) { for(int ir=0; ir < nxy; ir++) { @@ -754,7 +754,7 @@ void WF_atomic::atomicrandom(ModuleBase::ComplexMatrix& psi, for (int iw = iw_start ;iw < iw_end;iw++) { int startig = 0; - for(int ip = 0 ; ip < GlobalV::NPOL; ++ip) + for(int ip = 0 ; ip < PARAM.globalv.npol; ++ip) { for(int ig = 0 ; ig < npw ; ++ig) { diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp index 0f0b43f360..6984bfb33e 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp @@ -224,13 +224,13 @@ void Sto_Forces::cal_sto_force_nl(ModuleBase::matrix& forcenl, psi_in->fix_k(ik); stowf.shchi->fix_k(ik); //KS orbitals - int npmks = GlobalV::NPOL * nksbands; + int npmks = PARAM.globalv.npol * nksbands; zgemm_(&transa,&transb,&nkb,&npmks,&npw,&ModuleBase::ONE, GlobalC::ppcell.vkb.c,&npwx, psi_in->get_pointer(),&npwx, &ModuleBase::ZERO,becp.c,&nkb); //stochastic orbitals - int npmsto = GlobalV::NPOL * nstobands; + int npmsto = PARAM.globalv.npol * nstobands; zgemm_(&transa,&transb,&nkb,&npmsto,&npw,&ModuleBase::ONE, GlobalC::ppcell.vkb.c,&npwx, stowf.shchi->get_pointer(),&npwx, diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_hchi.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_hchi.cpp index 29a337f1ab..cc6b67015c 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_hchi.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_hchi.cpp @@ -1,5 +1,6 @@ #include "sto_hchi.h" +#include "module_parameter/parameter.h" #include "module_base/parallel_reduce.h" #include "module_base/timer.h" #include "module_base/tool_title.h" @@ -31,7 +32,7 @@ void Stochastic_hchi::hchi(complex* chig, complex* hchig, const const int current_spin = pkv->isk[ik]; const int npwx = this->wfcpw->npwk_max; const int npw = this->wfcpw->npwk[ik]; - const int npm = GlobalV::NPOL * m; + const int npm = PARAM.globalv.npol * m; const int inc = 1; const double tpiba2 = GlobalC::ucell.tpiba2; const int nrxx = this->wfcpw->nrxx; @@ -88,11 +89,11 @@ void Stochastic_hchi::hchi(complex* chig, complex* hchig, const if (GlobalC::ppcell.nkb > 0) { int nkb = GlobalC::ppcell.nkb; - complex* becp = new complex[nkb * GlobalV::NPOL * m]; + complex* becp = new complex[nkb * PARAM.globalv.npol * m]; char transc = 'C'; char transn = 'N'; char transt = 'T'; - if (m == 1 && GlobalV::NPOL == 1) + if (m == 1 && PARAM.globalv.npol == 1) { zgemv_(&transc, &npw, @@ -122,10 +123,10 @@ void Stochastic_hchi::hchi(complex* chig, complex* hchig, const becp, &nkb); } - Parallel_Reduce::reduce_pool(becp, nkb * GlobalV::NPOL * m); + Parallel_Reduce::reduce_pool(becp, nkb * PARAM.globalv.npol * m); - complex* Ps = new complex[nkb * GlobalV::NPOL * m]; - ModuleBase::GlobalFunc::ZEROS(Ps, GlobalV::NPOL * m * nkb); + complex* Ps = new complex[nkb * PARAM.globalv.npol * m]; + ModuleBase::GlobalFunc::ZEROS(Ps, PARAM.globalv.npol * m * nkb); int sum = 0; int iat = 0; @@ -153,7 +154,7 @@ void Stochastic_hchi::hchi(complex* chig, complex* hchig, const } // end na } // end nt - if (GlobalV::NPOL == 1 && m == 1) + if (PARAM.globalv.npol == 1 && m == 1) { zgemv_(&transn, &npw, diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp index ad5fa8feef..a0a00d8b9b 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp @@ -1,5 +1,6 @@ #include "sto_stress_pw.h" +#include "module_parameter/parameter.h" #include "module_base/timer.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" @@ -115,8 +116,9 @@ void Sto_Stress_PW::sto_stress_kin(ModuleBase::matrix& sigma, double twobysqrtpi = 2.0 / std::sqrt(ModuleBase::PI); double* kfac = new double[npwx]; int nksbands = psi_in->get_nbands(); - if (GlobalV::MY_STOGROUP != 0) + if (GlobalV::MY_STOGROUP != 0) { nksbands = 0; +} for (int ik = 0; ik < nks; ++ik) { @@ -236,8 +238,9 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma, int* nchip = stowf.nchip; const int npwx = wfc_basis->npwk_max; int nksbands = psi_in->get_nbands(); - if (GlobalV::MY_STOGROUP != 0) + if (GlobalV::MY_STOGROUP != 0) { nksbands = 0; +} // vkb1: |Beta(nkb,npw)> ModuleBase::ComplexMatrix vkb1(nkb, npwx); @@ -274,7 +277,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma, psi_in[0].fix_k(ik); stowf.shchi->fix_k(ik); // KS orbitals - int npmks = GlobalV::NPOL * nksbands; + int npmks = PARAM.globalv.npol * nksbands; zgemm_(&transa, &transb, &nkb, @@ -289,7 +292,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma, becp.c, &nkb); // stochastic orbitals - int npmsto = GlobalV::NPOL * nstobands; + int npmsto = PARAM.globalv.npol * nstobands; zgemm_(&transa, &transb, &nkb, @@ -363,10 +366,11 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma, { qvec = wfc_basis->getgpluskcar(ik, ig); double qm1; - if (qvec.norm2() > 1e-16) + if (qvec.norm2() > 1e-16) { qm1 = 1.0 / qvec.norm(); - else + } else { qm1 = 0; +} pdbecp_noevc[ig] -= 2.0 * pvkb[ig] * qvec0[ipol][0] * qvec0[jpol][0] * qm1 * this->ucell->tpiba; } // end ig } // end i @@ -403,10 +407,11 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma, for (int ib = 0; ib < nbandstot; ++ib) { double fac; - if (ib < nksbands) + if (ib < nksbands) { fac = wg(ik, ib); - else + } else { fac = p_kv->wk[ik]; +} int iat = 0; int sum = 0; for (int it = 0; it < this->ucell->ntype; ++it) diff --git a/source/module_hsolver/diago_elpa_native.cpp b/source/module_hsolver/diago_elpa_native.cpp index fb59f7024e..fa73a1d1e0 100644 --- a/source/module_hsolver/diago_elpa_native.cpp +++ b/source/module_hsolver/diago_elpa_native.cpp @@ -2,6 +2,7 @@ #include "module_base/blacs_connector.h" #include "module_base/global_variable.h" +#include "module_parameter/parameter.h" #include "module_base/lapack_connector.h" #include "module_base/scalapack_connector.h" #include "module_base/timer.h" @@ -94,7 +95,7 @@ void DiagoElpaNative::diag_pool(hamilt::MatrixBlock& h_mat, elpa_set(handle, "solver", ELPA_SOLVER_1STAGE, &success); #ifdef __CUDA - if (GlobalV::device_flag == "gpu") + if (PARAM.globalv.device_flag == "gpu") { elpa_set(handle, "nvidia-gpu", 1, &success); diff --git a/source/module_io/cal_r_overlap_R.cpp b/source/module_io/cal_r_overlap_R.cpp index cd07266507..3f4c35a6ba 100644 --- a/source/module_io/cal_r_overlap_R.cpp +++ b/source/module_io/cal_r_overlap_R.cpp @@ -309,12 +309,12 @@ void cal_r_overlap_R::out_rR(const int& istep) ic = this->ParaV->global2local_col(iw2); if (ic >= 0) { - int orb_index_row = iw1 / GlobalV::NPOL; - int orb_index_col = iw2 / GlobalV::NPOL; + int orb_index_row = iw1 / PARAM.globalv.npol; + int orb_index_col = iw2 / PARAM.globalv.npol; // The off-diagonal term in SOC calculaiton is zero, and the two diagonal terms are the same - int new_index = iw1 - GlobalV::NPOL * orb_index_row - + (iw2 - GlobalV::NPOL * orb_index_col) * GlobalV::NPOL; + int new_index = iw1 - PARAM.globalv.npol * orb_index_row + + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; if (new_index == 0 || new_index == 3) { @@ -583,12 +583,12 @@ void cal_r_overlap_R::out_rR_other(const int& istep, const std::setParaV->global2local_col(iw2); if (ic >= 0) { - int orb_index_row = iw1 / GlobalV::NPOL; - int orb_index_col = iw2 / GlobalV::NPOL; + int orb_index_row = iw1 / PARAM.globalv.npol; + int orb_index_col = iw2 / PARAM.globalv.npol; // The off-diagonal term in SOC calculaiton is zero, and the two diagonal terms are the same - int new_index = iw1 - GlobalV::NPOL * orb_index_row - + (iw2 - GlobalV::NPOL * orb_index_col) * GlobalV::NPOL; + int new_index = iw1 - PARAM.globalv.npol * orb_index_row + + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; if (new_index == 0 || new_index == 3) { diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 827ed1567d..abe21ba929 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -289,37 +289,6 @@ void Input_Conv::Convert() } } #endif - //-------------------------------------------- - // added by zhengdy-soc - //-------------------------------------------- - if (PARAM.inp.noncolin || PARAM.inp.lspinorb) - { - GlobalV::NSPIN = 4; - } - - if (GlobalV::NSPIN == 4) - { - GlobalV::NONCOLIN = PARAM.inp.noncolin; - // wavefunctions are spinors with 2 components - GlobalV::NPOL = 2; - // set the domag variable to make a spin-orbit calculation with zero - // magnetization - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = true; - GlobalV::LSPINORB = PARAM.inp.lspinorb; - if (PARAM.globalv.gamma_only_local) - { - ModuleBase::WARNING_QUIT("input_conv", - "nspin=4(soc or noncollinear-spin) does " - "not support gamma only calculation"); - } - } else { - GlobalV::LSPINORB = false; - GlobalV::NONCOLIN = false; - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; - GlobalV::NPOL = 1; - } //---------------------------------------------------------- // Yu Liu add 2022-05-18 @@ -488,7 +457,7 @@ void Input_Conv::Convert() } // In these case, inversion symmetry is also not allowed, symmetry should be // reset to -1 - if (GlobalV::LSPINORB) + if (PARAM.inp.lspinorb) { ModuleSymmetry::Symmetry::symm_flag = -1; } @@ -514,7 +483,6 @@ void Input_Conv::Convert() // wavefunction / charge / potential / (2/4) //---------------------------------------------------------- GlobalV::nelec = PARAM.inp.nelec; - GlobalV::out_pot = PARAM.inp.out_pot; #ifdef __LCAO diff --git a/source/module_io/read_set_globalv.cpp b/source/module_io/read_set_globalv.cpp index 1981cb9106..3760cc51fb 100644 --- a/source/module_io/read_set_globalv.cpp +++ b/source/module_io/read_set_globalv.cpp @@ -49,6 +49,48 @@ void ReadInput::set_globalv(Parameter& para) { para.sys.deepks_setorb = true; } + switch (para.input.nspin) + { + case 4: + if (para.input.noncolin) + { + para.sys.domag = true; + para.sys.domag_z = false; + } + else + { + para.sys.domag = false; + para.sys.domag_z = true; + } + para.sys.npol = 2; + break; + case 2: + case 1: + para.sys.domag = false; + para.sys.domag_z = false; + para.sys.npol = 1; + default: + break; + } + + if (para.input.device == "cpu") + { + para.sys.device_flag = "cpu"; + } + else if (para.input.device == "gpu") + { + if (para.input.basis_type == "lcao_in_pw") + { + GlobalV::ofs_warning << "The GPU currently does not support the basis type \"lcao_in_pw\"!" << std::endl; + para.sys.device_flag = "cpu"; + } + para.sys.device_flag = "gpu"; + } + else + { + GlobalV::ofs_warning << "Parameter \"device\" can only be set to \"cpu\" or \"gpu\"!" << std::endl; + ModuleBase::WARNING_QUIT("device", "Parameter \"device\" can only be set to \"cpu\" or \"gpu\"!"); + } } } @@ -74,6 +116,13 @@ void ReadInput::set_globalv_bcast() add_string_bcast(sys.global_matrix_dir); add_bool_bcast(sys.deepks_setorb); + + add_bool_bcast(sys.domag); + add_bool_bcast(sys.domag_z); + add_int_bcast(sys.npol); + + add_string_bcast(sys.device_flag); + add_bool_bcast(sys.double_grid); add_double_bcast(sys.uramping); } diff --git a/source/module_io/read_wfc_pw.cpp b/source/module_io/read_wfc_pw.cpp index 020f00f11e..3ed60c2d6b 100644 --- a/source/module_io/read_wfc_pw.cpp +++ b/source/module_io/read_wfc_pw.cpp @@ -1,5 +1,6 @@ #include "read_wfc_pw.h" +#include "module_parameter/parameter.h" #include "binstream.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" @@ -81,7 +82,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, ikstot = ik; #endif - npwtot *= GlobalV::NPOL; + npwtot *= PARAM.globalv.npol; @@ -200,7 +201,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, { glo_order[i] = -1; } - for (int i = 0; i < npwtot_in / GlobalV::NPOL; ++i) + for (int i = 0; i < npwtot_in / PARAM.globalv.npol; ++i) { int index = (miller[i].x * ny + miller[i].y) * nz + miller[i].z; glo_order[index] = i; @@ -255,7 +256,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, ip + GlobalV::NPROC_IN_POOL, POOL_WORLD, MPI_STATUS_IGNORE); - if (GlobalV::NPOL == 2) + if (PARAM.globalv.npol == 2) { MPI_Recv(&wfc(ib, npwk_max), pw_wfc->npwk[ik], @@ -280,7 +281,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, wfc_ip[i] = wfc_in[glo_order[ig_ip[i]]]; } MPI_Send(wfc_ip, size, MPI_DOUBLE_COMPLEX, ip, ip + GlobalV::NPROC_IN_POOL, POOL_WORLD); - if (GlobalV::NPOL == 2) + if (PARAM.globalv.npol == 2) { for (int i = 0; i < size; i++) { @@ -300,7 +301,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, { wfc(ib, i) = wfc_in[glo_order[l2g_pw[i]]]; } - if (GlobalV::NPOL == 2) + if (PARAM.globalv.npol == 2) { for (int i = 0; i < pw_wfc->npwk[ik]; ++i) { @@ -316,7 +317,7 @@ void ModuleIO::read_wfc_pw(const std::string& filename, { wfc(ib, i) = wfc_in[glo_order[l2g_pw[i]]]; } - if (GlobalV::NPOL == 2) + if (PARAM.globalv.npol == 2) { for (int i = 0; i < pw_wfc->npwk[ik]; ++i) { diff --git a/source/module_io/test/write_orb_info_test.cpp b/source/module_io/test/write_orb_info_test.cpp index 6d17f15749..0916f1444d 100644 --- a/source/module_io/test/write_orb_info_test.cpp +++ b/source/module_io/test/write_orb_info_test.cpp @@ -45,7 +45,7 @@ TEST(OrbInfo,WriteOrbInfo) ofs.open("running.log"); PARAM.sys.global_out_dir = "./"; PARAM.input.pseudo_rcut = 15.0; - GlobalV::LSPINORB = false; + PARAM.input.lspinorb = false; GlobalV::NSPIN = 1; PARAM.input.basis_type = "pw"; PARAM.input.dft_functional = "default"; diff --git a/source/module_io/test_serial/io_system_variable_test.cpp b/source/module_io/test_serial/io_system_variable_test.cpp index 7385246a77..b02395afc6 100644 --- a/source/module_io/test_serial/io_system_variable_test.cpp +++ b/source/module_io/test_serial/io_system_variable_test.cpp @@ -60,6 +60,22 @@ TEST_F(InputTest, Item_test) param.input.deepks_out_labels = true; readinput.set_globalv(param); EXPECT_EQ(param.sys.deepks_setorb, 1); + + param.input.nspin = 4; + param.input.noncolin = true; + readinput.set_globalv(param); + EXPECT_EQ(param.sys.domag, 1); + EXPECT_EQ(param.sys.domag_z, 0); + EXPECT_EQ(param.sys.npol, 2); + + param.input.nspin = 1; + param.input.lspinorb = true; + param.input.noncolin = false; + readinput.set_globalv(param); + EXPECT_EQ(param.sys.domag, 0); + EXPECT_EQ(param.sys.domag_z, 0); + EXPECT_EQ(param.sys.npol, 1); + } } \ No newline at end of file diff --git a/source/module_io/to_wannier90_lcao.cpp b/source/module_io/to_wannier90_lcao.cpp index a98f1c2ce3..1223c86503 100644 --- a/source/module_io/to_wannier90_lcao.cpp +++ b/source/module_io/to_wannier90_lcao.cpp @@ -688,10 +688,10 @@ void toWannier90_LCAO::construct_overlap_table_project() int row = this->ParaV->get_row_size(); int global_ir = 0; - for (int ir = 0; ir < row; ir += GlobalV::NPOL) + for (int ir = 0; ir < row; ir += PARAM.globalv.npol) { global_ir = ParaV->local2global_row(ir); - int orb_index_row = global_ir / GlobalV::NPOL; + int orb_index_row = global_ir / PARAM.globalv.npol; for (int wannier_index = 0; wannier_index < num_wannier; wannier_index++) { @@ -767,7 +767,7 @@ void toWannier90_LCAO::cal_orbA_overlap_R() for (int ir = 0; ir < row; ir++) { global_ir = ParaV->local2global_row(ir); - int orb_index_row = global_ir / GlobalV::NPOL; + int orb_index_row = global_ir / PARAM.globalv.npol; int it1 = iw2it[orb_index_row]; int ia1 = iw2ia[orb_index_row]; diff --git a/source/module_io/to_wannier90_lcao_in_pw.cpp b/source/module_io/to_wannier90_lcao_in_pw.cpp index 17c14d8b04..9cc6bb9acf 100644 --- a/source/module_io/to_wannier90_lcao_in_pw.cpp +++ b/source/module_io/to_wannier90_lcao_in_pw.cpp @@ -1,5 +1,6 @@ #include "to_wannier90_lcao_in_pw.h" +#include "module_parameter/parameter.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_base/math_integral.h" #include "module_base/math_polyint.h" @@ -113,7 +114,7 @@ psi::Psi>* toWannier90_LCAO_IN_PW::get_unk_from_lcao( { // init int npwx = wfcpw->npwk_max; - psi::Psi> *unk_inLcao = new psi::Psi>(num_kpts, num_bands, npwx*GlobalV::NPOL, kv.ngk.data()); + psi::Psi> *unk_inLcao = new psi::Psi>(num_kpts, num_bands, npwx*PARAM.globalv.npol, kv.ngk.data()); unk_inLcao->zero_out(); // Orbital projection to plane wave @@ -122,7 +123,7 @@ psi::Psi>* toWannier90_LCAO_IN_PW::get_unk_from_lcao( for (int ik = 0; ik < num_kpts; ik++) { int npw = kv.ngk[ik]; - ModuleBase::ComplexMatrix orbital_in_G(GlobalV::NLOCAL, npwx*GlobalV::NPOL); + ModuleBase::ComplexMatrix orbital_in_G(GlobalV::NLOCAL, npwx*PARAM.globalv.npol); // Wavefunc_in_pw::produce_local_basis_in_pw(ik, wfcpw, sf, orbital_in_G, table_local); //produce_local_basis_in_pw(ik, wfcpw, sf, orbital_in_G, table_local); nao_G_expansion(ik, wfcpw, orbital_in_G); @@ -162,7 +163,7 @@ psi::Psi>* toWannier90_LCAO_IN_PW::get_unk_from_lcao( { for (int ib = 0; ib < num_bands; ib++) { - // for (int ig = 0; ig < npwx*GlobalV::NPOL; ig++) + // for (int ig = 0; ig < npwx*PARAM.globalv.npol; ig++) // { // for (int iw = 0; iw < GlobalV::NLOCAL; iw++) // { @@ -219,7 +220,7 @@ void toWannier90_LCAO_IN_PW::nao_G_expansion( if(psig.expired()) { ModuleBase::WARNING_QUIT("toWannier90_LCAO_IN_PW::nao_G_expansion", "psig is expired"); } int nbands = GlobalV::NLOCAL; - int nbasis = npwx*GlobalV::NPOL; + int nbasis = npwx*PARAM.globalv.npol; for (int ib = 0; ib < nbands; ib++) { for (int ig = 0; ig < nbasis; ig++) diff --git a/source/module_io/unk_overlap_pw.cpp b/source/module_io/unk_overlap_pw.cpp index c1477e8b89..58dfdd9f97 100644 --- a/source/module_io/unk_overlap_pw.cpp +++ b/source/module_io/unk_overlap_pw.cpp @@ -1,5 +1,6 @@ #include "unk_overlap_pw.h" +#include "module_parameter/parameter.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" unkOverlap_pw::unkOverlap_pw() @@ -135,12 +136,12 @@ std::complex unkOverlap_pw::unkdotp_soc_G(const ModulePW::PW_Basis_K* wf std::complex result(0.0,0.0); const int number_pw = wfcpw->npw; - std::complex* unk_L = new std::complex[number_pw * GlobalV::NPOL]; - std::complex* unk_R = new std::complex[number_pw * GlobalV::NPOL]; - ModuleBase::GlobalFunc::ZEROS(unk_L,number_pw*GlobalV::NPOL); - ModuleBase::GlobalFunc::ZEROS(unk_R,number_pw*GlobalV::NPOL); + std::complex* unk_L = new std::complex[number_pw * PARAM.globalv.npol]; + std::complex* unk_R = new std::complex[number_pw * PARAM.globalv.npol]; + ModuleBase::GlobalFunc::ZEROS(unk_L,number_pw*PARAM.globalv.npol); + ModuleBase::GlobalFunc::ZEROS(unk_R,number_pw*PARAM.globalv.npol); - for(int i = 0; i < GlobalV::NPOL; i++) + for(int i = 0; i < PARAM.globalv.npol; i++) { for (int igl = 0; igl < evc->get_ngk(ik_L); igl++) { @@ -153,7 +154,7 @@ std::complex unkOverlap_pw::unkdotp_soc_G(const ModulePW::PW_Basis_K* wf } } - for (int iG = 0; iG < number_pw*GlobalV::NPOL; iG++) + for (int iG = 0; iG < number_pw*PARAM.globalv.npol; iG++) { result = result + conj(unk_L[iG]) * unk_R[iG]; @@ -220,7 +221,7 @@ std::complex unkOverlap_pw::unkdotp_soc_G0(const ModulePW::PW_Basis* rho wfcpw->real2recip(psi_up, psi_up, ik_L); wfcpw->real2recip(psi_down, psi_down, ik_L); - for (int i = 0; i < GlobalV::NPOL; i++) + for (int i = 0; i < PARAM.globalv.npol; i++) { for(int ig = 0; ig < evc->get_ngk(ik_R); ig++) { diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index d409d66a13..c44c3109cc 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -1,5 +1,6 @@ #include "pot_hxc_lrtd.h" #include "module_elecstate/potentials/pot_base.h" +#include "module_parameter/parameter.h" #include "module_elecstate/potentials/H_Hartree_pw.h" #include "module_base/timer.h" #include "module_hamilt_general/module_xc/xc_functional.h" @@ -14,7 +15,7 @@ namespace LR { this->rho_basis_ = rho_basis_in; this->nrxx = chg_gs->nrxx; - this->nspin = (GlobalV::NSPIN == 1 || (GlobalV::NSPIN == 4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) ? 1 : 2; + this->nspin = (GlobalV::NSPIN == 1 || (GlobalV::NSPIN == 4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) ? 1 : 2; this->pot_hartree = LR_Util::make_unique(this->rho_basis_); std::set local_xc = { "lda", "pbe", "hse" }; diff --git a/source/module_parameter/system_parameter.h b/source/module_parameter/system_parameter.h index 72061a9794..f4a1cac33b 100644 --- a/source/module_parameter/system_parameter.h +++ b/source/module_parameter/system_parameter.h @@ -38,8 +38,11 @@ struct System_para std::string global_matrix_dir = ""; ///< global matrix directory bool deepks_setorb = false; ///< true if "deepks" is set - - + int npol = 1; ///< number of polarization + bool domag = false; /// 1 : calculate the magnetism with x, y, z component + bool domag_z = false; /// 1 : constrain the magnetism to z axis + + std::string device_flag = "cpu"; ///< device flag, "cpu" or "gpu" bool double_grid = false; ///< true if "ndx,ndy,ndz" is larger than "nx,ny,nz" double uramping = -10.0 / 13.6; /// U-Ramping method (Ry) std::vector hubbard_u = {}; ///< Hubbard Coulomb interaction parameter U (Ry) diff --git a/source/module_psi/psi.cpp b/source/module_psi/psi.cpp index c026b237ac..a164da46a8 100644 --- a/source/module_psi/psi.cpp +++ b/source/module_psi/psi.cpp @@ -1,5 +1,6 @@ #include "psi.h" +#include "module_parameter/parameter.h" #include "module_base/global_variable.h" #include "module_base/tool_quit.h" @@ -31,7 +32,7 @@ Range::Range(const bool k_first_in, const size_t index_1_in, const size_t range_ template Psi::Psi() { - this->npol = GlobalV::NPOL; + this->npol = PARAM.globalv.npol; this->device = base_device::get_device_type(this->ctx); } @@ -43,7 +44,7 @@ template Psi::~Psi() template Psi::Psi(const int* ngk_in) { this->ngk = ngk_in; - this->npol = GlobalV::NPOL; + this->npol = PARAM.globalv.npol; this->device = base_device::get_device_type(this->ctx); } @@ -53,7 +54,7 @@ template Psi::Psi(const int nk_in, cons this->ngk = ngk_in; this->current_b = 0; this->current_k = 0; - this->npol = GlobalV::NPOL; + this->npol = PARAM.globalv.npol; this->device = base_device::get_device_type(this->ctx); this->resize(nk_in, nbd_in, nbs_in); // Currently only GPU's implementation is supported for device recording! @@ -70,7 +71,7 @@ template Psi::Psi(T* psi_pointer, const this->ngk = ngk_in; this->current_b = 0; this->current_k = 0; - this->npol = GlobalV::NPOL; + this->npol = PARAM.globalv.npol; this->device = base_device::get_device_type(this->ctx); this->nk = nk_in; this->nbands = nbd_in; diff --git a/source/module_psi/psi_initializer.cpp b/source/module_psi/psi_initializer.cpp index d214a13da3..6773b11241 100644 --- a/source/module_psi/psi_initializer.cpp +++ b/source/module_psi/psi_initializer.cpp @@ -53,20 +53,20 @@ psi::Psi>* psi_initializer::allocate(bool only_p { /* EVERY ZETA FOR (2l+1) ORBS */ /* - non-rotate basis, nbands_local*=2 for GlobalV::NPOL = 2 is enough + non-rotate basis, nbands_local*=2 for PARAM.globalv.npol = 2 is enough */ - //nbands_local += this->p_ucell_->atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL; + //nbands_local += this->p_ucell_->atoms[it].l_nchi[l]*(2*l+1) * PARAM.globalv.npol; /* rotate basis, nbands_local*=4 for p, d, f,... orbitals, and nbands_local*=2 for s orbitals risky when NSPIN = 4, problematic psi value, needed to be checked */ if(l == 0) { - nbands_local += this->p_ucell_->atoms[it].l_nchi[l] * GlobalV::NPOL; + nbands_local += this->p_ucell_->atoms[it].l_nchi[l] * PARAM.globalv.npol; } else { - nbands_local += this->p_ucell_->atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL; + nbands_local += this->p_ucell_->atoms[it].l_nchi[l]*(2*l+1) * PARAM.globalv.npol; } } } @@ -84,7 +84,7 @@ psi::Psi>* psi_initializer::allocate(bool only_p } } int nkpts_actual = (PARAM.inp.calculation == "nscf" && this->mem_saver_ == 1)? 1 : this->pw_wfc_->nks; - int nbasis_actual = this->pw_wfc_->npwk_max * GlobalV::NPOL; + int nbasis_actual = this->pw_wfc_->npwk_max * PARAM.globalv.npol; psi::Psi>* psi_out = nullptr; if(!only_psig) { @@ -97,7 +97,7 @@ psi::Psi>* psi_initializer::allocate(bool only_p */ const size_t memory_cost_psi = nkpts_actual* - GlobalV::NBANDS * this->pw_wfc_->npwk_max * GlobalV::NPOL* + GlobalV::NBANDS * this->pw_wfc_->npwk_max * PARAM.globalv.npol* sizeof(std::complex); std::cout << " MEMORY FOR PSI PER PROCESSOR (MB) : " << double(memory_cost_psi)/1024.0/1024.0 << std::endl; ModuleBase::Memory::record("Psi_PW", memory_cost_psi); @@ -108,7 +108,7 @@ psi::Psi>* psi_initializer::allocate(bool only_p this->pw_wfc_->npwk); const size_t memory_cost_psig = nkpts_actual* - nbands_actual * this->pw_wfc_->npwk_max * GlobalV::NPOL* + nbands_actual * this->pw_wfc_->npwk_max * PARAM.globalv.npol* sizeof(T); std::cout << " MEMORY FOR AUXILLARY PSI PER PROCESSOR (MB) : " << double(memory_cost_psig)/1024.0/1024.0 << std::endl; @@ -122,7 +122,7 @@ psi::Psi>* psi_initializer::allocate(bool only_p << "nbands_complem = " << this->nbands_complem_ << "\n" << "nbasis_actual = " << nbasis_actual << "\n" << "npwk_max = " << this->pw_wfc_->npwk_max << "\n" - << "npol = " << GlobalV::NPOL << "\n"; + << "npol = " << PARAM.globalv.npol << "\n"; ModuleBase::Memory::record("psigPW", memory_cost_psig); ModuleBase::timer::tick("psi_initializer", "allocate"); return psi_out; @@ -149,9 +149,9 @@ void psi_initializer::random_t(T* psi, const int iw_start, const int for (int iw = iw_start; iw < iw_end; iw++) { // get the starting memory address of iw band - T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * GlobalV::NPOL]); + T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * PARAM.globalv.npol]); int startig = 0; - for(int ipol = 0 ; ipol < GlobalV::NPOL ; ++ipol) + for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) { for(int ir=0; ir < nxy; ir++) @@ -191,7 +191,7 @@ void psi_initializer::random_t(T* psi, const int iw_start, const int #endif for (int iw = iw_start ;iw < iw_end; iw++) { - T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * GlobalV::NPOL]); + T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * PARAM.globalv.npol]); for (int ig = 0; ig < ng; ig++) { const double rr = std::rand()/double(RAND_MAX); //qianrui add RAND_MAX @@ -199,7 +199,7 @@ void psi_initializer::random_t(T* psi, const int iw_start, const int const double gk2 = this->pw_wfc_->getgk2(ik,ig); psi_slice[ig] = this->template cast_to_T(std::complex(rr*cos(arg)/(gk2 + 1.0), rr*sin(arg)/(gk2 + 1.0))); } - if(GlobalV::NPOL==2) + if(PARAM.globalv.npol==2) { for (int ig = this->pw_wfc_->npwk_max; ig < this->pw_wfc_->npwk_max + ng; ig++) { diff --git a/source/module_psi/psi_initializer_atomic.cpp b/source/module_psi/psi_initializer_atomic.cpp index b89fe3e62f..a10963f5d2 100644 --- a/source/module_psi/psi_initializer_atomic.cpp +++ b/source/module_psi/psi_initializer_atomic.cpp @@ -202,7 +202,7 @@ void psi_initializer_atomic::proj_ao_onkG(int ik) Soc soc; soc.rot_ylm(l + 1); const double j = this->p_ucell_->atoms[it].ncpp.jchi[ipswfc]; /* NOT NONCOLINEAR CASE, rotation matrix become identity */ - if (!(GlobalV::DOMAG||GlobalV::DOMAG_Z)) + if (!(PARAM.globalv.domag||PARAM.globalv.domag_z)) { double cg_coeffs[2]; for(int m = -l-1; m < l+1; m++) diff --git a/source/module_psi/psi_initializer_nao.cpp b/source/module_psi/psi_initializer_nao.cpp index 8a5e96b74a..bd0b2dfe37 100644 --- a/source/module_psi/psi_initializer_nao.cpp +++ b/source/module_psi/psi_initializer_nao.cpp @@ -467,7 +467,7 @@ void psi_initializer_nao::proj_ao_onkG(int ik) ibasis += 2 * L + 1; } } // end for is_N - } // end if GlobalV::NONCOLIN + } // end if PARAM.inp.noncolin else { // LSDA and nomagnet case /* DOES NOT DISTINGUISH m QUANTUM NUMBER FOR CHI */ diff --git a/source/module_psi/test/psi_initializer_unit_test.cpp b/source/module_psi/test/psi_initializer_unit_test.cpp index df170684a5..8dceaa2a8f 100644 --- a/source/module_psi/test/psi_initializer_unit_test.cpp +++ b/source/module_psi/test/psi_initializer_unit_test.cpp @@ -119,12 +119,12 @@ class PsiIntializerUnitTest : public ::testing::Test { GlobalV::NSPIN = 1; PARAM.input.orbital_dir = "./support/"; PARAM.input.pseudo_dir = "./support/"; - GlobalV::NPOL = 1; + PARAM.sys.npol = 1; PARAM.input.calculation = "scf"; PARAM.input.init_wfc = "random"; GlobalV::KS_SOLVER = "cg"; - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; // lattice this->p_ucell->a1 = {10.0, 0.0, 0.0}; this->p_ucell->a2 = {0.0, 10.0, 0.0}; @@ -525,7 +525,7 @@ TEST_F(PsiIntializerUnitTest, CalPsigAtomic) { TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc) { PARAM.input.init_wfc = "atomic"; GlobalV::NSPIN = 4; - GlobalV::NPOL = 2; + PARAM.sys.npol = 2; this->p_ucell->atoms[0].ncpp.has_so = false; this->p_ucell->natomwfc *= 2; this->psi_init = new psi_initializer_atomic, base_device::DEVICE_CPU>(); @@ -549,7 +549,7 @@ TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc) { this->psi_init->proj_ao_onkG(0); EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); GlobalV::NSPIN = 1; - GlobalV::NPOL = 1; + PARAM.sys.npol = 1; this->p_ucell->atoms[0].ncpp.has_so = false; this->p_ucell->natomwfc /= 2; delete psi; @@ -558,7 +558,7 @@ TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc) { TEST_F(PsiIntializerUnitTest, CalPsigAtomicSocHasSo) { PARAM.input.init_wfc = "atomic"; GlobalV::NSPIN = 4; - GlobalV::NPOL = 2; + PARAM.sys.npol = 2; this->p_ucell->atoms[0].ncpp.has_so = true; this->p_ucell->natomwfc *= 2; this->psi_init = new psi_initializer_atomic, base_device::DEVICE_CPU>(); @@ -582,7 +582,7 @@ TEST_F(PsiIntializerUnitTest, CalPsigAtomicSocHasSo) { this->psi_init->proj_ao_onkG(0); EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); GlobalV::NSPIN = 1; - GlobalV::NPOL = 1; + PARAM.sys.npol = 1; this->p_ucell->atoms[0].ncpp.has_so = false; this->p_ucell->natomwfc /= 2; delete psi; @@ -666,10 +666,10 @@ TEST_F(PsiIntializerUnitTest, CalPsigNaoRandom) { TEST_F(PsiIntializerUnitTest, CalPsigNaoSoc) { PARAM.input.init_wfc = "nao"; GlobalV::NSPIN = 4; - GlobalV::NPOL = 2; + PARAM.sys.npol = 2; this->p_ucell->atoms[0].ncpp.has_so = false; - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; this->psi_init = new psi_initializer_nao, base_device::DEVICE_CPU>(); #ifdef __MPI this->psi_init->initialize(this->p_sf, @@ -696,10 +696,10 @@ TEST_F(PsiIntializerUnitTest, CalPsigNaoSoc) { TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSo) { PARAM.input.init_wfc = "nao"; GlobalV::NSPIN = 4; - GlobalV::NPOL = 2; + PARAM.sys.npol = 2; this->p_ucell->atoms[0].ncpp.has_so = true; - GlobalV::DOMAG = false; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = false; + PARAM.sys.domag_z = false; this->psi_init = new psi_initializer_nao, base_device::DEVICE_CPU>(); #ifdef __MPI this->psi_init->initialize(this->p_sf, @@ -726,10 +726,10 @@ TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSo) { TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSoDOMAG) { PARAM.input.init_wfc = "nao"; GlobalV::NSPIN = 4; - GlobalV::NPOL = 2; + PARAM.sys.npol = 2; this->p_ucell->atoms[0].ncpp.has_so = true; - GlobalV::DOMAG = true; - GlobalV::DOMAG_Z = false; + PARAM.sys.domag = true; + PARAM.sys.domag_z = false; this->psi_init = new psi_initializer_nao, base_device::DEVICE_CPU>(); #ifdef __MPI this->psi_init->initialize(this->p_sf, From 865773f3177d089b48ec59a83017cd23681d25fe Mon Sep 17 00:00:00 2001 From: LUNASEA <33978601+maki49@users.noreply.github.com> Date: Wed, 18 Sep 2024 20:44:35 +0800 Subject: [PATCH 2/8] enlarge threshold (#5117) --- tests/integrate/281_NO_KP_HSE_symmetry/threshold | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/integrate/281_NO_KP_HSE_symmetry/threshold b/tests/integrate/281_NO_KP_HSE_symmetry/threshold index 03285b49d0..b213a90bd4 100644 --- a/tests/integrate/281_NO_KP_HSE_symmetry/threshold +++ b/tests/integrate/281_NO_KP_HSE_symmetry/threshold @@ -1,4 +1,4 @@ -threshold 0.00001 -force_threshold 0.00001 -stress_threshold 0.00001 -fatal_threshold 0.00001 \ No newline at end of file +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.001 +fatal_threshold 0.001 \ No newline at end of file From 026a4b578023caabf4ebb0c25bddc1c37e9c3b22 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 18 Sep 2024 20:44:57 +0800 Subject: [PATCH 3/8] Build(deps): Bump pre-commit-ci/lite-action from 1.0.2 to 1.0.3 (#5110) Bumps [pre-commit-ci/lite-action](https://github.com/pre-commit-ci/lite-action) from 1.0.2 to 1.0.3. - [Release notes](https://github.com/pre-commit-ci/lite-action/releases) - [Commits](https://github.com/pre-commit-ci/lite-action/compare/v1.0.2...v1.0.3) --- updated-dependencies: - dependency-name: pre-commit-ci/lite-action dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index de425874df..2f9996e070 100644 --- a/.github/workflows/test.yml +++ b/.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.2 + - uses: pre-commit-ci/lite-action@v1.0.3 - name: Build run: | From a9e89e4d42725ccffb3d6668c012b7ac60ebd2a0 Mon Sep 17 00:00:00 2001 From: Haozhi Han Date: Wed, 18 Sep 2024 20:49:55 +0800 Subject: [PATCH 4/8] Refactor: refactor paw code in HSolverLIP (#5126) * refactor paw code in HSolverLIP * fix build bug * fix build bug --- source/module_hsolver/hsolver_lcaopw.cpp | 455 ++++++++++++----------- source/module_hsolver/hsolver_lcaopw.h | 78 ++-- 2 files changed, 273 insertions(+), 260 deletions(-) diff --git a/source/module_hsolver/hsolver_lcaopw.cpp b/source/module_hsolver/hsolver_lcaopw.cpp index 4291a9d713..c1261ff695 100644 --- a/source/module_hsolver/hsolver_lcaopw.cpp +++ b/source/module_hsolver/hsolver_lcaopw.cpp @@ -1,10 +1,5 @@ #include "hsolver_lcaopw.h" -#include "module_parameter/parameter.h" -#include "diago_bpcg.h" -#include "diago_cg.h" -#include "diago_dav_subspace.h" -#include "diago_david.h" #include "module_base/global_variable.h" #include "module_base/parallel_global.h" // for MPI #include "module_base/timer.h" @@ -15,8 +10,8 @@ #include "module_hamilt_pw/hamilt_pwdft/wavefunc.h" #include "module_hsolver/diagh.h" #include "module_hsolver/diago_iter_assist.h" +#include "module_parameter/parameter.h" -#include #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" #endif @@ -24,263 +19,281 @@ #ifdef __EXX #include "module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h" #endif + namespace hsolver { - template - HSolverLIP::HSolverLIP(ModulePW::PW_Basis_K* wfc_basis_in) +#ifdef USE_PAW +template +void HSolverLIP::paw_func_in_kloop(const int ik) +{ + if (PARAM.inp.use_paw) { - this->wfc_basis = wfc_basis_in; - } + const int npw = this->wfc_basis->npwk[ik]; + ModuleBase::Vector3* _gk = new ModuleBase::Vector3[npw]; + for (int ig = 0; ig < npw; ig++) + { + _gk[ig] = this->wfc_basis->getgpluskcar(ik, ig); + } + + std::vector kpt(3, 0); + kpt[0] = this->wfc_basis->kvec_c[ik].x; + kpt[1] = this->wfc_basis->kvec_c[ik].y; + kpt[2] = this->wfc_basis->kvec_c[ik].z; + + double** kpg; + double** gcar; + kpg = new double*[npw]; + gcar = new double*[npw]; + for (int ipw = 0; ipw < npw; ipw++) + { + kpg[ipw] = new double[3]; + kpg[ipw][0] = _gk[ipw].x; + kpg[ipw][1] = _gk[ipw].y; + kpg[ipw][2] = _gk[ipw].z; + + gcar[ipw] = new double[3]; + gcar[ipw][0] = this->wfc_basis->getgcar(ik, ipw).x; + gcar[ipw][1] = this->wfc_basis->getgcar(ik, ipw).y; + gcar[ipw][2] = this->wfc_basis->getgcar(ik, ipw).z; + } + + GlobalC::paw_cell.set_paw_k(npw, + wfc_basis->npwk_max, + kpt.data(), + this->wfc_basis->get_ig2ix(ik).data(), + this->wfc_basis->get_ig2iy(ik).data(), + this->wfc_basis->get_ig2iz(ik).data(), + (const double**)kpg, + GlobalC::ucell.tpiba, + (const double**)gcar); + + std::vector().swap(kpt); + for (int ipw = 0; ipw < npw; ipw++) + { + delete[] kpg[ipw]; + delete[] gcar[ipw]; + } + delete[] kpg; + delete[] gcar; + GlobalC::paw_cell.get_vkb(); - /* - lcao_in_pw - */ - template - void HSolverLIP::solve(hamilt::Hamilt* pHamilt, // ESolver_KS_PW::p_hamilt - psi::Psi& psi, // ESolver_KS_PW::kspw_psi - elecstate::ElecState* pes, // ESolver_KS_PW::pes - psi::Psi& transform, - const bool skip_charge) + GlobalC::paw_cell.set_currentk(ik); + } +} + +template +void HSolverLIP::paw_func_after_kloop(psi::Psi& psi, elecstate::ElecState* pes) +{ + if (PARAM.inp.use_paw) { - ModuleBase::TITLE("HSolverLIP", "solve"); - ModuleBase::timer::tick("HSolverLIP", "solve"); - std::vector eigenvalues(pes->ekb.nr * pes->ekb.nc, 0); + if (typeid(Real) != typeid(double)) + { + ModuleBase::WARNING_QUIT("HSolverLIP::solve", "PAW is only supported for double precision!"); + } + + GlobalC::paw_cell.reset_rhoij(); for (int ik = 0; ik < this->wfc_basis->nks; ++ik) { - /// update H(k) for each k point - pHamilt->updateHk(ik); -#ifdef USE_PAW - if (PARAM.inp.use_paw) + const int npw = this->wfc_basis->npwk[ik]; + ModuleBase::Vector3* _gk = new ModuleBase::Vector3[npw]; + for (int ig = 0; ig < npw; ig++) { - const int npw = this->wfc_basis->npwk[ik]; - ModuleBase::Vector3* _gk = new ModuleBase::Vector3[npw]; - for (int ig = 0; ig < npw; ig++) - { - _gk[ig] = this->wfc_basis->getgpluskcar(ik, ig); - } + _gk[ig] = this->wfc_basis->getgpluskcar(ik, ig); + } - std::vector kpt(3, 0); - kpt[0] = this->wfc_basis->kvec_c[ik].x; - kpt[1] = this->wfc_basis->kvec_c[ik].y; - kpt[2] = this->wfc_basis->kvec_c[ik].z; + std::vector kpt(3, 0); + kpt[0] = this->wfc_basis->kvec_c[ik].x; + kpt[1] = this->wfc_basis->kvec_c[ik].y; + kpt[2] = this->wfc_basis->kvec_c[ik].z; - double** kpg; - double** gcar; - kpg = new double* [npw]; - gcar = new double* [npw]; - for (int ipw = 0; ipw < npw; ipw++) - { - kpg[ipw] = new double[3]; - kpg[ipw][0] = _gk[ipw].x; - kpg[ipw][1] = _gk[ipw].y; - kpg[ipw][2] = _gk[ipw].z; - - gcar[ipw] = new double[3]; - gcar[ipw][0] = this->wfc_basis->getgcar(ik, ipw).x; - gcar[ipw][1] = this->wfc_basis->getgcar(ik, ipw).y; - gcar[ipw][2] = this->wfc_basis->getgcar(ik, ipw).z; - } + double** kpg; + double** gcar; + kpg = new double*[npw]; + gcar = new double*[npw]; + for (int ipw = 0; ipw < npw; ipw++) + { + kpg[ipw] = new double[3]; + kpg[ipw][0] = _gk[ipw].x; + kpg[ipw][1] = _gk[ipw].y; + kpg[ipw][2] = _gk[ipw].z; + + gcar[ipw] = new double[3]; + gcar[ipw][0] = this->wfc_basis->getgcar(ik, ipw).x; + gcar[ipw][1] = this->wfc_basis->getgcar(ik, ipw).y; + gcar[ipw][2] = this->wfc_basis->getgcar(ik, ipw).z; + } - GlobalC::paw_cell.set_paw_k(npw, - wfc_basis->npwk_max, - kpt.data(), - this->wfc_basis->get_ig2ix(ik).data(), - this->wfc_basis->get_ig2iy(ik).data(), - this->wfc_basis->get_ig2iz(ik).data(), - (const double**)kpg, - GlobalC::ucell.tpiba, - (const double**)gcar); - - std::vector().swap(kpt); - for (int ipw = 0; ipw < npw; ipw++) - { - delete[] kpg[ipw]; - delete[] gcar[ipw]; - } - delete[] kpg; - delete[] gcar; + GlobalC::paw_cell.set_paw_k(npw, + wfc_basis->npwk_max, + kpt.data(), + this->wfc_basis->get_ig2ix(ik).data(), + this->wfc_basis->get_ig2iy(ik).data(), + this->wfc_basis->get_ig2iz(ik).data(), + (const double**)kpg, + GlobalC::ucell.tpiba, + (const double**)gcar); + + std::vector().swap(kpt); + for (int ipw = 0; ipw < npw; ipw++) + { + delete[] kpg[ipw]; + delete[] gcar[ipw]; + } + delete[] kpg; + delete[] gcar; - GlobalC::paw_cell.get_vkb(); + GlobalC::paw_cell.get_vkb(); - GlobalC::paw_cell.set_currentk(ik); - } -#endif psi.fix_k(ik); - transform.fix_k(ik); - /// solve eigenvector and eigenvalue for H(k) - - // hsolver::DiagoIterAssist::diagH_subspace( - // pHamilt, // interface to hamilt - // transform, // transform matrix between lcao and pw - // psi, // psi in pw basis - // eigenvalues.data() + ik * pes->ekb.nc, // eigenvalues - // psi.get_nbands() // number of the lowest energies bands - // ); -#ifdef __EXX - auto& exx_lip = dynamic_cast*>(pHamilt)->exx_lip; - auto add_exx_to_subspace_hamilt = [&ik, &exx_lip](T* hcc, const int naos) -> void - { - if (GlobalC::exx_info.info_global.cal_exx) { - for (int n = 0; n < naos; ++n) { - for (int m = 0; m < naos; ++m) { - hcc[n * naos + m] += (T)GlobalC::exx_info.info_global.hybrid_alpha * - exx_lip.get_exx_matrix()[ik][m][n]; -} -} -} - }; - auto set_exxlip_lcaowfc = [&ik, &exx_lip](const T* const vcc, const int naos, const int nbands) -> void - { - if (GlobalC::exx_info.info_global.cal_exx) { - exx_lip.set_hvec(ik, vcc, naos, nbands); -} - }; -#endif - hsolver::DiagoIterAssist::diagH_subspace_init( - pHamilt, // interface to hamilt - transform.get_pointer(), // transform matrix between lcao and pw - transform.get_nbands(), - transform.get_nbasis(), - psi, // psi in pw basis - eigenvalues.data() + ik * pes->ekb.nc // eigenvalues -#ifdef __EXX - , add_exx_to_subspace_hamilt - , set_exxlip_lcaowfc -#endif - ); - - if (skip_charge) + GlobalC::paw_cell.set_currentk(ik); + int nbands = psi.get_nbands(); + for (int ib = 0; ib < nbands; ib++) { - GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik - << " is: " << DiagoIterAssist::avg_iter - << " ; where current threshold is: " << DiagoIterAssist::PW_DIAG_THR - << " . " << std::endl; - DiagoIterAssist::avg_iter = 0.0; + GlobalC::paw_cell.accumulate_rhoij(reinterpret_cast*>(psi.get_pointer(ib)), + pes->wg(ik, ib)); } - /// calculate the contribution of Psi for charge density rho } - castmem_2d_2h_op()(cpu_ctx, cpu_ctx, pes->ekb.c, eigenvalues.data(), pes->ekb.nr * pes->ekb.nc); - if (skip_charge) - { - ModuleBase::timer::tick("HSolverLIP", "solve"); - return; - } - reinterpret_cast*>(pes)->psiToRho(psi); + std::vector> rhoijp; + std::vector> rhoijselect; + std::vector nrhoijsel; -#ifdef USE_PAW - if (PARAM.inp.use_paw) +#ifdef __MPI + if (GlobalV::RANK_IN_POOL == 0) { - if (typeid(Real) != typeid(double)) - { - ModuleBase::WARNING_QUIT("HSolverLIP::solve", "PAW is only supported for double precision!"); - } + GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - GlobalC::paw_cell.reset_rhoij(); - for (int ik = 0; ik < this->wfc_basis->nks; ++ik) + for (int iat = 0; iat < GlobalC::ucell.nat; iat++) { - const int npw = this->wfc_basis->npwk[ik]; - ModuleBase::Vector3* _gk = new ModuleBase::Vector3[npw]; - for (int ig = 0; ig < npw; ig++) - { - _gk[ig] = this->wfc_basis->getgpluskcar(ik, ig); - } + GlobalC::paw_cell.set_rhoij(iat, + nrhoijsel[iat], + rhoijselect[iat].size(), + rhoijselect[iat].data(), + rhoijp[iat].data()); + } + } +#else + GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - std::vector kpt(3, 0); - kpt[0] = this->wfc_basis->kvec_c[ik].x; - kpt[1] = this->wfc_basis->kvec_c[ik].y; - kpt[2] = this->wfc_basis->kvec_c[ik].z; + for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + { + GlobalC::paw_cell.set_rhoij(iat, + nrhoijsel[iat], + rhoijselect[iat].size(), + rhoijselect[iat].data(), + rhoijp[iat].data()); + } - double** kpg; - double** gcar; - kpg = new double* [npw]; - gcar = new double* [npw]; - for (int ipw = 0; ipw < npw; ipw++) - { - kpg[ipw] = new double[3]; - kpg[ipw][0] = _gk[ipw].x; - kpg[ipw][1] = _gk[ipw].y; - kpg[ipw][2] = _gk[ipw].z; - - gcar[ipw] = new double[3]; - gcar[ipw][0] = this->wfc_basis->getgcar(ik, ipw).x; - gcar[ipw][1] = this->wfc_basis->getgcar(ik, ipw).y; - gcar[ipw][2] = this->wfc_basis->getgcar(ik, ipw).z; - } +#endif + double* nhatgr; + GlobalC::paw_cell.get_nhat(pes->charge->nhat, nhatgr); + } +} +#endif - GlobalC::paw_cell.set_paw_k(npw, - wfc_basis->npwk_max, - kpt.data(), - this->wfc_basis->get_ig2ix(ik).data(), - this->wfc_basis->get_ig2iy(ik).data(), - this->wfc_basis->get_ig2iz(ik).data(), - (const double**)kpg, - GlobalC::ucell.tpiba, - (const double**)gcar); - - std::vector().swap(kpt); - for (int ipw = 0; ipw < npw; ipw++) - { - delete[] kpg[ipw]; - delete[] gcar[ipw]; - } - delete[] kpg; - delete[] gcar; +template +HSolverLIP::HSolverLIP(ModulePW::PW_Basis_K* wfc_basis_in) +{ + this->wfc_basis = wfc_basis_in; +} - GlobalC::paw_cell.get_vkb(); +/* + lcao_in_pw +*/ +template +void HSolverLIP::solve(hamilt::Hamilt* pHamilt, // ESolver_KS_PW::p_hamilt + psi::Psi& psi, // ESolver_KS_PW::kspw_psi + elecstate::ElecState* pes, // ESolver_KS_PW::pes + psi::Psi& transform, + const bool skip_charge) +{ + ModuleBase::TITLE("HSolverLIP", "solve"); + ModuleBase::timer::tick("HSolverLIP", "solve"); + std::vector eigenvalues(pes->ekb.nr * pes->ekb.nc, 0); + for (int ik = 0; ik < this->wfc_basis->nks; ++ik) + { + /// update H(k) for each k point + pHamilt->updateHk(ik); - psi.fix_k(ik); - GlobalC::paw_cell.set_currentk(ik); - int nbands = psi.get_nbands(); - for (int ib = 0; ib < nbands; ib++) - { - GlobalC::paw_cell.accumulate_rhoij(reinterpret_cast*>(psi.get_pointer(ib)), - pes->wg(ik, ib)); - } - } +#ifdef USE_PAW + this->paw_func_in_kloop(ik); +#endif - std::vector> rhoijp; - std::vector> rhoijselect; - std::vector nrhoijsel; + psi.fix_k(ik); + transform.fix_k(ik); -#ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) +#ifdef __EXX + auto& exx_lip = dynamic_cast*>(pHamilt)->exx_lip; + auto add_exx_to_subspace_hamilt = [&ik, &exx_lip](T* hcc, const int naos) -> void { + if (GlobalC::exx_info.info_global.cal_exx) { - GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int n = 0; n < naos; ++n) { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); + for (int m = 0; m < naos; ++m) + { + hcc[n * naos + m] + += (T)GlobalC::exx_info.info_global.hybrid_alpha * exx_lip.get_exx_matrix()[ik][m][n]; + } } } -#else - GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + }; + auto set_exxlip_lcaowfc = [&ik, &exx_lip](const T* const vcc, const int naos, const int nbands) -> void { + if (GlobalC::exx_info.info_global.cal_exx) { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); + exx_lip.set_hvec(ik, vcc, naos, nbands); } - + }; #endif - double* nhatgr; - GlobalC::paw_cell.get_nhat(pes->charge->nhat, nhatgr); - } + /// solve eigenvector and eigenvalue for H(k) + hsolver::DiagoIterAssist::diagH_subspace_init( + pHamilt, // interface to hamilt + transform.get_pointer(), // transform matrix between lcao and pw + transform.get_nbands(), + transform.get_nbasis(), + psi, // psi in pw basis + eigenvalues.data() + ik * pes->ekb.nc // eigenvalues +#ifdef __EXX + , + add_exx_to_subspace_hamilt, + set_exxlip_lcaowfc #endif + ); + + if (skip_charge) + { + GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik + << " is: " << DiagoIterAssist::avg_iter + << " ; where current threshold is: " << DiagoIterAssist::PW_DIAG_THR << " . " + << std::endl; + DiagoIterAssist::avg_iter = 0.0; + } + /// calculate the contribution of Psi for charge density rho + } + base_device::memory::cast_memory_op()( + cpu_ctx, + cpu_ctx, + pes->ekb.c, + eigenvalues.data(), + pes->ekb.nr * pes->ekb.nc); + + if (skip_charge) + { ModuleBase::timer::tick("HSolverLIP", "solve"); return; } + reinterpret_cast*>(pes)->psiToRho(psi); + +#ifdef USE_PAW + this->paw_func_after_kloop(psi, pes); +#endif + + ModuleBase::timer::tick("HSolverLIP", "solve"); + return; +} - template class HSolverLIP>; - template class HSolverLIP>; +template class HSolverLIP>; +template class HSolverLIP>; } // namespace hsolver diff --git a/source/module_hsolver/hsolver_lcaopw.h b/source/module_hsolver/hsolver_lcaopw.h index a606721d67..1694e6c62b 100644 --- a/source/module_hsolver/hsolver_lcaopw.h +++ b/source/module_hsolver/hsolver_lcaopw.h @@ -4,45 +4,45 @@ #include "hsolver.h" #include "module_base/macros.h" #include "module_base/module_device/types.h" -namespace hsolver { - - // LCAO-in-PW does not support GPU now. - template - class HSolverLIP - { - private: - // Note GetTypeReal::type will - // return T if T is real type(float, double), - // otherwise return the real type of T(complex, complex) - using Real = typename GetTypeReal::type; - - public: - - HSolverLIP(ModulePW::PW_Basis_K* wfc_basis_in); - - /// @brief solve function for lcao_in_pw - /// @param pHamilt interface to hamilt - /// @param psi reference to psi - /// @param pes interface to elecstate - /// @param transform transformation matrix between lcao and pw - /// @param skip_charge - void solve(hamilt::Hamilt* pHamilt, - psi::Psi& psi, - elecstate::ElecState* pes, - psi::Psi& transform, - const bool skip_charge); - - protected: - - ModulePW::PW_Basis_K* wfc_basis = nullptr; - - std::vector eigenvalues; - using castmem_2d_2h_op - = base_device::memory::cast_memory_op; - - private: - std::string method = "none"; - }; +namespace hsolver +{ + +// LCAO-in-PW does not support GPU now. +template +class HSolverLIP +{ + private: + // Note GetTypeReal::type will + // return T if T is real type(float, double), + // otherwise return the real type of T(complex, complex) + using Real = typename GetTypeReal::type; + + public: + HSolverLIP(ModulePW::PW_Basis_K* wfc_basis_in); + + /// @brief solve function for lcao_in_pw + /// @param pHamilt interface to hamilt + /// @param psi reference to psi + /// @param pes interface to elecstate + /// @param transform transformation matrix between lcao and pw + /// @param skip_charge + void solve(hamilt::Hamilt* pHamilt, + psi::Psi& psi, + elecstate::ElecState* pes, + psi::Psi& transform, + const bool skip_charge); + + private: + ModulePW::PW_Basis_K* wfc_basis = nullptr; + + std::vector eigenvalues; + +#ifdef USE_PAW + void paw_func_in_kloop(const int ik); + + void paw_func_after_kloop(psi::Psi& psi, elecstate::ElecState* pes); +#endif +}; } // namespace hsolver From 43809b8786149ab54562589cf3b36008906716b0 Mon Sep 17 00:00:00 2001 From: Haozhi Han Date: Wed, 18 Sep 2024 20:50:43 +0800 Subject: [PATCH 5/8] remove my_conj func from hsolver_lcao (#5128) --- source/module_hsolver/hsolver_lcao.cpp | 4 ++-- source/module_hsolver/hsolver_lcao.h | 12 ------------ 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/source/module_hsolver/hsolver_lcao.cpp b/source/module_hsolver/hsolver_lcao.cpp index 6d7692ccf3..66818d6eb7 100644 --- a/source/module_hsolver/hsolver_lcao.cpp +++ b/source/module_hsolver/hsolver_lcao.cpp @@ -208,8 +208,8 @@ void HSolverLCAO::hamiltSolvePsiK(hamilt::Hamilt* hm, psi::Psi& { for (int j = i; j < h_mat.col; j++) { - h_mat.p[h_mat.row * j + i] = hsolver::my_conj(h_mat.p[h_mat.row * i + j]); - s_mat.p[s_mat.row * j + i] = hsolver::my_conj(s_mat.p[s_mat.row * i + j]); + h_mat.p[h_mat.row * j + i] = hsolver::get_conj(h_mat.p[h_mat.row * i + j]); + s_mat.p[s_mat.row * j + i] = hsolver::get_conj(s_mat.p[s_mat.row * i + j]); } } diff --git a/source/module_hsolver/hsolver_lcao.h b/source/module_hsolver/hsolver_lcao.h index 9d27f50c19..8dbb3e52f8 100644 --- a/source/module_hsolver/hsolver_lcao.h +++ b/source/module_hsolver/hsolver_lcao.h @@ -36,18 +36,6 @@ class HSolverLCAO std::string method = "none"; }; -template -inline T my_conj(T value) -{ - return value; -} - -template <> -inline std::complex my_conj(std::complex value) -{ - return std::conj(value); -} - } // namespace hsolver #endif \ No newline at end of file From 7d4391773b50786a809e5b2143fd21137adc47c1 Mon Sep 17 00:00:00 2001 From: kirk0830 <67682086+kirk0830@users.noreply.github.com> Date: Thu, 19 Sep 2024 10:32:51 +0800 Subject: [PATCH 6/8] Refactor&Fix: optimize the memory usage by psi_initializer (#5120) * change some non-const input to const * fix a severe bug that random number cannot fill initialized bands when init_wfc nao * resolve conflicts by hand --- .../module_hamilt_pw/hamilt_pwdft/wfinit.cpp | 30 ++--- source/module_hsolver/diago_dav_subspace.cpp | 2 +- source/module_io/print_info.cpp | 35 ++---- source/module_psi/psi_initializer.cpp | 118 ++++++++++-------- source/module_psi/psi_initializer.h | 6 +- source/module_psi/psi_initializer_atomic.cpp | 37 ++++-- source/module_psi/psi_initializer_atomic.h | 2 +- .../psi_initializer_atomic_random.cpp | 11 +- .../psi_initializer_atomic_random.h | 2 +- source/module_psi/psi_initializer_nao.cpp | 13 +- source/module_psi/psi_initializer_nao.h | 2 +- .../module_psi/psi_initializer_nao_random.cpp | 11 +- .../module_psi/psi_initializer_nao_random.h | 2 +- source/module_psi/psi_initializer_random.cpp | 5 +- source/module_psi/psi_initializer_random.h | 2 +- 15 files changed, 148 insertions(+), 130 deletions(-) diff --git a/source/module_hamilt_pw/hamilt_pwdft/wfinit.cpp b/source/module_hamilt_pw/hamilt_pwdft/wfinit.cpp index 04332f7795..a29e646f08 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/wfinit.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/wfinit.cpp @@ -136,18 +136,7 @@ void WFInit::allocate_psi(Psi>*& psi, template void WFInit::make_table(const int nks, Structure_Factor* p_sf) { - if (this->use_psiinitializer) // new initialization method, used in KSDFT and LCAO_IN_PW calculation - { - // re-tabulate because GlobalV::DQ may change due to the change of atomic positions and cell parameters - // for nao, we recalculate the overlap matrix between flz and jlq - // for atomic, we recalculate the overlap matrix between pswfc and jlq - // for psig is not read-only, its value will be overwritten in initialize_psi(), dont need delete and - // reallocate - if ((this->init_wfc.substr(0, 3) == "nao") || (this->init_wfc.substr(0, 6) == "atomic")) - { - this->psi_init->tabulate(); - } - } + if (this->use_psiinitializer) {} // do not need to do anything because the interpolate table is unchanged else // old initialization method, used in EXX calculation { this->p_wf->init_after_vc(nks); // reallocate wanf2, the planewave expansion of lcao @@ -161,11 +150,13 @@ void WFInit::initialize_psi(Psi>* psi, hamilt::Hamilt* p_hamilt, std::ofstream& ofs_running) { - if (!this->use_psiinitializer) { - return; -} + if (!this->use_psiinitializer) { return; } ModuleBase::timer::tick("WFInit", "initialize_psi"); + // if psig is not allocated before, allocate it if (!this->psi_init->psig_use_count()) { this->psi_init->allocate(/*psig_only=*/true); } + + // loop over kpoints, make it possible to only allocate memory for psig at the only one kpt + // like (1, nbands, npwx), in which npwx is the maximal npw of all kpoints for (int ik = 0; ik < this->pw_wfc->nks; ik++) { //! Fix the wavefunction to initialize at given kpoint @@ -192,7 +183,12 @@ void WFInit::initialize_psi(Psi>* psi, //! to use psig, we need to lock it to get a shared pointer version, //! then switch kpoint of psig to the given one auto psig_ = psig.lock(); - psig_->fix_k(ik); + // CHANGE LOG: if not lcaoinpw, the psig will only be used in psi-initialization + // so we can only allocate memory for one kpoint with the maximal number of pw + // over all kpoints, then the memory space will be always enough. Then for each + // kpoint, the psig is calculated in an overwrite manner. + const int ik_psig = (psig_->get_nk() == 1) ? 0 : ik; + psig_->fix_k(ik_psig); std::vector::type> etatom(psig_->get_nbands(), 0.0); @@ -200,7 +196,7 @@ void WFInit::initialize_psi(Psi>* psi, // either by matrix-multiplication or by copying-discarding if (this->psi_init->method() != "random") { - // lcao_in_pw and pw share the same esolver. In the future, we will have different esolver + // lcaoinpw and pw share the same esolver. In the future, we will have different esolver if (((this->ks_solver == "cg") || (this->ks_solver == "lapack")) && (this->basis_type == "pw")) { // the following function is only run serially, to be improved diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp index f6d99d8447..99e0065137 100644 --- a/source/module_hsolver/diago_dav_subspace.cpp +++ b/source/module_hsolver/diago_dav_subspace.cpp @@ -816,7 +816,7 @@ int Diago_DavSubspace::diag(const HPsiFunc& hpsi_func, do { - printf("enter diag... is_subspace = %d, ntry = %d\n", this->is_subspace, ntry); + //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); diff --git a/source/module_io/print_info.cpp b/source/module_io/print_info.cpp index 13960c9f91..2d15c07701 100644 --- a/source/module_io/print_info.cpp +++ b/source/module_io/print_info.cpp @@ -74,10 +74,9 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) << std::setw(16) << "KPOINTS" << std::setw(12) << "PROCESSORS"; - if(PARAM.inp.basis_type=="lcao" || PARAM.inp.basis_type=="lcao_in_pw" || (PARAM.inp.basis_type=="pw" && PARAM.inp.init_wfc.substr(0, 3) == "nao")) - { - std::cout << std::setw(12) << "NBASE"; - } + const bool orbinfo = (PARAM.inp.basis_type=="lcao" || PARAM.inp.basis_type=="lcao_in_pw" + || (PARAM.inp.basis_type=="pw" && PARAM.inp.init_wfc.substr(0, 3) == "nao")); + if (orbinfo) { std::cout << std::setw(12) << "NBASE"; } std::cout << std::endl; std::cout << " " << std::setw(8) << GlobalV::NSPIN; @@ -92,11 +91,7 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) } std::cout << std::setw(12) << GlobalV::NPROC; - - if(PARAM.inp.basis_type=="lcao" || PARAM.inp.basis_type=="lcao_in_pw" || (PARAM.inp.basis_type=="pw" && PARAM.inp.init_wfc.substr(0, 3) == "nao")) - { - std::cout << std::setw(12) << GlobalV::NLOCAL; - } + if (orbinfo) { std::cout << std::setw(12) << GlobalV::NLOCAL; } std::cout << std::endl; @@ -104,15 +99,15 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) std::cout << " ---------------------------------------------------------" << std::endl; - if(PARAM.inp.basis_type=="lcao") + if(PARAM.inp.basis_type == "lcao") { std::cout << " Use Systematically Improvable Atomic bases" << std::endl; } - else if(PARAM.inp.basis_type=="lcao_in_pw") + else if(PARAM.inp.basis_type == "lcao_in_pw") { std::cout << " Expand Atomic bases into plane waves" << std::endl; } - else if(PARAM.inp.basis_type=="pw") + else if(PARAM.inp.basis_type == "pw") { std::cout << " Use plane wave basis" << std::endl; } @@ -126,7 +121,7 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) std::cout << " " << std::setw(8) << "ELEMENT"; - if(PARAM.inp.basis_type=="lcao" || PARAM.inp.basis_type=="lcao_in_pw") + if (orbinfo) { std::cout << std::setw(16) << "ORBITALS"; std::cout << std::setw(12) << "NBASE"; @@ -137,29 +132,21 @@ void Print_Info::setup_parameters(UnitCell &ucell, K_Vectors &kv) std::cout << std::endl; - + const std::string spectrum = "spdfghi"; for(int it=0; it -psi::Psi>* psi_initializer::allocate(bool only_psig) +psi::Psi>* psi_initializer::allocate(const bool only_psig) { ModuleBase::timer::tick("psi_initializer", "allocate"); /* @@ -83,40 +88,48 @@ psi::Psi>* psi_initializer::allocate(bool only_p } } } - int nkpts_actual = (PARAM.inp.calculation == "nscf" && this->mem_saver_ == 1)? 1 : this->pw_wfc_->nks; - int nbasis_actual = this->pw_wfc_->npwk_max * PARAM.globalv.npol; + + const int nks_psi = (PARAM.inp.calculation == "nscf" && this->mem_saver_ == 1)? 1 : this->pw_wfc_->nks; + const int nbasis_actual = this->pw_wfc_->npwk_max * PARAM.globalv.npol; psi::Psi>* psi_out = nullptr; if(!only_psig) { - psi_out = new psi::Psi>(nkpts_actual, + psi_out = new psi::Psi>(nks_psi, GlobalV::NBANDS, // because no matter what, the wavefunction finally needed has GlobalV::NBANDS bands nbasis_actual, this->pw_wfc_->npwk); - /* - WARNING: this will cause DIRECT MEMORY LEAK, psi is not properly deallocated - */ - const size_t memory_cost_psi = - nkpts_actual* - GlobalV::NBANDS * this->pw_wfc_->npwk_max * PARAM.globalv.npol* + double memory_cost_psi = nks_psi * GlobalV::NBANDS * this->pw_wfc_->npwk_max * PARAM.globalv.npol* sizeof(std::complex); - std::cout << " MEMORY FOR PSI PER PROCESSOR (MB) : " << double(memory_cost_psi)/1024.0/1024.0 << std::endl; +#ifdef __MPI + // get the correct memory cost for psi by all-reduce sum + Parallel_Reduce::reduce_all(memory_cost_psi); +#endif + // std::cout << " MEMORY FOR PSI PER PROCESSOR (MB) : " << double(memory_cost_psi)/1024.0/1024.0 << std::endl; ModuleBase::Memory::record("Psi_PW", memory_cost_psi); } - this->psig_ = std::make_shared>(nkpts_actual, + // for memory saving, the psig can always only hold one k-point data. But for lcao_in_pw, the psig + // is actcually a transformation matrix. During the SCF, the projection might be quite time- + // consuming. + const int nks_psig = (this->mem_saver_ == 1 && PARAM.inp.basis_type != "lcao_in_pw")? 1 : nks_psi; + this->psig_ = std::make_shared>(nks_psig, nbands_actual, nbasis_actual, this->pw_wfc_->npwk); - const size_t memory_cost_psig = - nkpts_actual* - nbands_actual * this->pw_wfc_->npwk_max * PARAM.globalv.npol* - sizeof(T); - std::cout << " MEMORY FOR AUXILLARY PSI PER PROCESSOR (MB) : " << double(memory_cost_psig)/1024.0/1024.0 << std::endl; + + double memory_cost_psig = + nks_psig * nbands_actual * this->pw_wfc_->npwk_max * PARAM.globalv.npol * sizeof(T); +#ifdef __MPI + // get the correct memory cost for psig by all-reduce sum + Parallel_Reduce::reduce_all(memory_cost_psig); +#endif + // std::cout << " MEMORY FOR AUXILLARY PSI PER PROCESSOR (MB) : " << double(memory_cost_psig)/1024.0/1024.0 << std::endl; GlobalV::ofs_running << "Allocate memory for psi and psig done.\n" << "Print detailed information of dimension of psi and psig:\n" - << "psi: (" << nkpts_actual << ", " << GlobalV::NBANDS << ", " << nbasis_actual << ")\n" - << "psig: (" << nkpts_actual << ", " << nbands_actual << ", " << nbasis_actual << ")\n" - << "nkpts_actual = " << nkpts_actual << "\n" + << "psi: (" << nks_psi << ", " << GlobalV::NBANDS << ", " << nbasis_actual << ")\n" + << "psig: (" << nks_psig << ", " << nbands_actual << ", " << nbasis_actual << ")\n" + << "nks (psi) = " << nks_psi << "\n" + << "nks (psig) = " << nks_psig << "\n" << "GlobalV::NBANDS = " << GlobalV::NBANDS << "\n" << "nbands_actual = " << nbands_actual << "\n" << "nbands_complem = " << this->nbands_complem_ << "\n" @@ -134,10 +147,15 @@ void psi_initializer::random_t(T* psi, const int iw_start, const int ModuleBase::timer::tick("psi_initializer", "random_t"); assert(iw_start >= 0); const int ng = this->pw_wfc_->npwk[ik]; -#ifdef __MPI + + // if random seed is specified, then based on this seed to generate random wavefunction if (this->random_seed_ > 0) // qianrui add 2021-8-13 { +#ifdef __MPI // if compile with MPI, then let the seed include the kpoint information srand(unsigned(this->random_seed_ + this->p_parakpts_->startk_pool[GlobalV::MY_POOL] + ik)); +#else // otherwise, it is the run in serial, without the Parallel_Kpoints class + srand(unsigned(this->random_seed_ + ik)); +#endif const int nxy = this->pw_wfc_->fftnxy; const int nz = this->pw_wfc_->nz; const int nstnz = this->pw_wfc_->nst*nz; @@ -153,66 +171,68 @@ void psi_initializer::random_t(T* psi, const int iw_start, const int int startig = 0; for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) { - - for(int ir=0; ir < nxy; ir++) + // loop over all fft (x,y), but actually loop over all sticks + for(int ir = 0; ir < nxy; ir++) { - if(this->pw_wfc_->fftixy2ip[ir] < 0) { continue; -} - if(GlobalV::RANK_IN_POOL==0) + // if the stick is not on present processor, then skip + if(this->pw_wfc_->fftixy2ip[ir] < 0) { continue; } + // otherwise + if(GlobalV::RANK_IN_POOL == 0) { - for(int iz=0; izpw_wfc_->getigl2isz(ik,ig)]; - const double arg= ModuleBase::TWO_PI * tmparg[this->pw_wfc_->getigl2isz(ik,ig)]; - const double gk2 = this->pw_wfc_->getgk2(ik,ig); + // get the correct value of "rr" and "arg" by indexing map "getigl2isz" + const double rr = tmprr[this->pw_wfc_->getigl2isz(ik, ig)]; + const double arg = ModuleBase::TWO_PI * tmparg[this->pw_wfc_->getigl2isz(ik,ig)]; + // get the |G+k|^2 + const double gk2 = this->pw_wfc_->getgk2(ik, ig); + // initialize the wavefunction value with rr/(gk2 + 1.0) * exp(i*arg) psi_slice[ig+startig] = this->template cast_to_T(std::complex(rr*cos(arg)/(gk2 + 1.0), rr*sin(arg)/(gk2 + 1.0))); } - startig += this->pw_wfc_->npwk_max; + startig += this->pw_wfc_->npwk_max; // move to the next polarization } } } - else + else // if random seed is not specified, then generate random wavefunction directly { -#else // !__MPI - if (this->random_seed_ > 0) // qianrui add 2021-8-13 - { - srand(unsigned(this->random_seed_ + ik)); - } -#endif for (int iw = iw_start ;iw < iw_end; iw++) { - T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * PARAM.globalv.npol]); + T* psi_slice = &(psi[iw * this->pw_wfc_->npwk_max * PARAM.globalv.npol]); // get the memory to write directly. For nspin 4, nbasis*2 for (int ig = 0; ig < ng; ig++) { const double rr = std::rand()/double(RAND_MAX); //qianrui add RAND_MAX - const double arg= ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); - const double gk2 = this->pw_wfc_->getgk2(ik,ig); + const double arg = ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); + const double gk2 = this->pw_wfc_->getgk2(ik, ig); psi_slice[ig] = this->template cast_to_T(std::complex(rr*cos(arg)/(gk2 + 1.0), rr*sin(arg)/(gk2 + 1.0))); } - if(PARAM.globalv.npol==2) + if(PARAM.globalv.npol==2) // additionally for nspin 4... { for (int ig = this->pw_wfc_->npwk_max; ig < this->pw_wfc_->npwk_max + ng; ig++) { const double rr = std::rand()/double(RAND_MAX); - const double arg= ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); - const double gk2 = this->pw_wfc_->getgk2(ik,ig-this->pw_wfc_->npwk_max); + const double arg = ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); + const double gk2 = this->pw_wfc_->getgk2(ik, ig - this->pw_wfc_->npwk_max); psi_slice[ig] = this->template cast_to_T(std::complex(rr*cos(arg)/(gk2 + 1.0), rr*sin(arg)/(gk2 + 1.0))); } } } -#ifdef __MPI } -#endif ModuleBase::timer::tick("psi_initializer_random", "random_t"); } diff --git a/source/module_psi/psi_initializer.h b/source/module_psi/psi_initializer.h index e481fb1dd8..059a8b258b 100644 --- a/source/module_psi/psi_initializer.h +++ b/source/module_psi/psi_initializer.h @@ -89,7 +89,7 @@ class psi_initializer pseudopot_cell_vnl* = nullptr) = 0; //< nonlocal pseudopotential #endif /// @brief CENTRAL FUNCTION: allocate memory for psi - psi::Psi>* allocate(bool only_psig = false); //< if only allocate memory for psig + psi::Psi>* allocate(const bool only_psig = false); //< if only allocate memory for psig void random_t(T* psi, //< [out] psi const int iw_start, //< iw_start, starting band index @@ -105,7 +105,7 @@ class psi_initializer /// @brief CENTRAL FUNCTION: calculate the interpolate table virtual void tabulate() = 0; /// @brief CENTRAL FUNCTION: calculate projection of atomic radial function onto planewave basis BASED ON THE OVERLAP TABLE - virtual void proj_ao_onkG(int ik) = 0; + virtual void proj_ao_onkG(const int ik) = 0; // getter and setter UnitCell* p_ucell() const { return this->p_ucell_; } @@ -175,7 +175,7 @@ class psi_initializer // avoid memory leak std::shared_ptr> psig_; private: - int mem_saver_ = 0; + int mem_saver_ = 1; std::string method_ = "none"; int nbands_complem_ = 0; double random_mix_ = 0; diff --git a/source/module_psi/psi_initializer_atomic.cpp b/source/module_psi/psi_initializer_atomic.cpp index a10963f5d2..7c19153c7a 100644 --- a/source/module_psi/psi_initializer_atomic.cpp +++ b/source/module_psi/psi_initializer_atomic.cpp @@ -55,8 +55,11 @@ void psi_initializer_atomic::initialize(Structure_Factor* sf, const int& rank) { ModuleBase::timer::tick("psi_initializer_atomic", "initialize"); - if(p_pspot_nl == nullptr) ModuleBase::WARNING_QUIT("psi_initializer_atomic::initialize_only_once", - "pseudopot_cell_vnl object cannot be mullptr for atomic, quit."); + if(p_pspot_nl == nullptr) + { + ModuleBase::WARNING_QUIT("psi_initializer_atomic::initialize", + "pseudopot_cell_vnl object cannot be nullptr for atomic, quit."); + } // import this->sf_ = sf; this->pw_wfc_ = pw_wfc; @@ -66,6 +69,10 @@ void psi_initializer_atomic::initialize(Structure_Factor* sf, this->random_seed_ = random_seed; // allocate this->allocate_table(); + // then for generate random number to fill in the wavefunction + this->ixy2is_.clear(); + this->ixy2is_.resize(this->pw_wfc_->fftnxy); + this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); ModuleBase::timer::tick("psi_initializer_atomic", "initialize_only_once"); } #else @@ -77,8 +84,11 @@ void psi_initializer_atomic::initialize(Structure_Factor* sf, pseudopot_cell_vnl* p_pspot_nl) { ModuleBase::timer::tick("psi_initializer_atomic", "initialize"); - if(p_pspot_nl == nullptr) ModuleBase::WARNING_QUIT("psi_initializer_atomic::initialize_only_once", - "pseudopot_cell_vnl object cannot be mullptr for atomic, quit."); + if(p_pspot_nl == nullptr) + { + ModuleBase::WARNING_QUIT("psi_initializer_atomic::initialize", + "pseudopot_cell_vnl object cannot be nullptr for atomic, quit."); + } // import this->sf_ = sf; this->pw_wfc_ = pw_wfc; @@ -87,6 +97,10 @@ void psi_initializer_atomic::initialize(Structure_Factor* sf, this->random_seed_ = random_seed; // allocate this->allocate_table(); + // then for generate random number to fill in the wavefunction + this->ixy2is_.clear(); + this->ixy2is_.resize(this->pw_wfc_->fftnxy); + this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); ModuleBase::timer::tick("psi_initializer_atomic", "initialize_only_once"); } #endif @@ -122,7 +136,7 @@ void psi_initializer_atomic::tabulate() { int n_rgrid = (PARAM.inp.pseudo_mesh)?atom->ncpp.mesh:atom->ncpp.msh; std::vector pswfcr(n_rgrid); - for (int ir=0; irncpp.chi(ic, ir); + for (int ir=0; irncpp.chi(ic, ir); } normalize(n_rgrid, pswfcr, atom->ncpp.rab.data()); if (atom->ncpp.oc[ic] >= 0.0) // reasonable occupation number, but is it always true? { @@ -141,17 +155,18 @@ void psi_initializer_atomic::tabulate() std::complex phase_factor(double arg, int mode) { - if(mode == 1) return std::complex(cos(arg),0); - else if (mode == -1) return std::complex(0, sin(arg)); - else if (mode == 0) return std::complex(cos(arg), sin(arg)); - else return std::complex(1,0); + if(mode == 1) { return std::complex(cos(arg),0); } + else if (mode == -1) { return std::complex(0, sin(arg)); } + else if (mode == 0) { return std::complex(cos(arg), sin(arg)); } + else { return std::complex(1,0); } } template -void psi_initializer_atomic::proj_ao_onkG(int ik) +void psi_initializer_atomic::proj_ao_onkG(const int ik) { ModuleBase::timer::tick("psi_initializer_atomic", "proj_ao_onkG"); - this->psig_->fix_k(ik); + const int ik_psig = (this->psig_->get_nk() == 1) ? 0 : ik; + this->psig_->fix_k(ik_psig); //this->print_status(psi); const int npw = this->pw_wfc_->npwk[ik]; int lmax = this->p_ucell_->lmax_ppwf; diff --git a/source/module_psi/psi_initializer_atomic.h b/source/module_psi/psi_initializer_atomic.h index 3eb3596649..0bb13b604b 100644 --- a/source/module_psi/psi_initializer_atomic.h +++ b/source/module_psi/psi_initializer_atomic.h @@ -34,7 +34,7 @@ class psi_initializer_atomic : public psi_initializer #endif virtual void allocate_table() override; virtual void tabulate() override; - virtual void proj_ao_onkG(int ik) override; + virtual void proj_ao_onkG(const int ik) override; // additional getter std::vector pseudopot_files() const { return pseudopot_files_; } diff --git a/source/module_psi/psi_initializer_atomic_random.cpp b/source/module_psi/psi_initializer_atomic_random.cpp index 5a8b71ca9f..b502e95d72 100644 --- a/source/module_psi/psi_initializer_atomic_random.cpp +++ b/source/module_psi/psi_initializer_atomic_random.cpp @@ -11,9 +11,6 @@ void psi_initializer_atomic_random::initialize(Structure_Factor* sf, const int& rank) { psi_initializer_atomic::initialize(sf, pw_wfc, p_ucell, p_parakpts, random_seed, p_pspot_nl, rank); - this->ixy2is_.clear(); - this->ixy2is_.resize(this->pw_wfc_->fftnxy); - this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); } #else template @@ -24,17 +21,15 @@ void psi_initializer_atomic_random::initialize(Structure_Factor* sf, pseudopot_cell_vnl* p_pspot_nl) { psi_initializer_atomic::initialize(sf, pw_wfc, p_ucell, random_seed, p_pspot_nl); - this->ixy2is_.clear(); - this->ixy2is_.resize(this->pw_wfc_->fftnxy); - this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); } #endif template -void psi_initializer_atomic_random::proj_ao_onkG(int ik) +void psi_initializer_atomic_random::proj_ao_onkG(const int ik) { double rm = this->random_mix(); - this->psig_->fix_k(ik); + const int ik_psig = (this->psig_->get_nk() == 1) ? 0 : ik; + this->psig_->fix_k(ik_psig); psi_initializer_atomic::proj_ao_onkG(ik); psi::Psi psi_random(1, this->psig_->get_nbands(), this->psig_->get_nbasis(), nullptr); psi_random.fix_k(0); diff --git a/source/module_psi/psi_initializer_atomic_random.h b/source/module_psi/psi_initializer_atomic_random.h index d884c28fd7..88b2bcfb93 100644 --- a/source/module_psi/psi_initializer_atomic_random.h +++ b/source/module_psi/psi_initializer_atomic_random.h @@ -34,7 +34,7 @@ class psi_initializer_atomic_random : public psi_initializer_atomic pseudopot_cell_vnl* = nullptr) override;//< nonlocal pseudopotential #endif - virtual void proj_ao_onkG(int ik) override; + virtual void proj_ao_onkG(const int ik) override; virtual void tabulate() override {psi_initializer_atomic::tabulate();}; private: diff --git a/source/module_psi/psi_initializer_nao.cpp b/source/module_psi/psi_initializer_nao.cpp index bd0b2dfe37..c831b59015 100644 --- a/source/module_psi/psi_initializer_nao.cpp +++ b/source/module_psi/psi_initializer_nao.cpp @@ -300,6 +300,10 @@ void psi_initializer_nao::initialize(Structure_Factor* sf, this->read_external_orbs(this->p_ucell_->orbital_fn, rank); // this->cal_ovlp_flzjlq(); //because GlobalV::NQX will change during vcrelax, so it should be called in both init // and init_after_vc + // then for generate random number to fill in the wavefunction + this->ixy2is_.clear(); + this->ixy2is_.resize(this->pw_wfc_->fftnxy); + this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); ModuleBase::timer::tick("psi_initializer_nao", "initialize_mpi"); } #else @@ -322,6 +326,10 @@ void psi_initializer_nao::initialize(Structure_Factor* sf, this->read_external_orbs(this->p_ucell_->orbital_fn, 0); // this->cal_ovlp_flzjlq(); //because GlobalV::NQX will change during vcrelax, so it should be called in both init // and init_after_vc + // then for generate random number to fill in the wavefunction + this->ixy2is_.clear(); + this->ixy2is_.resize(this->pw_wfc_->fftnxy); + this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); ModuleBase::timer::tick("psi_initializer_nao", "initialize_serial"); } #endif @@ -363,11 +371,12 @@ void psi_initializer_nao::tabulate() } template -void psi_initializer_nao::proj_ao_onkG(int ik) +void psi_initializer_nao::proj_ao_onkG(const int ik) { ModuleBase::timer::tick("psi_initializer_nao", "initialize"); assert(ik >= 0); - this->psig_->fix_k(ik); + const int ik_psig = (this->psig_->get_nk() == 1) ? 0 : ik; + this->psig_->fix_k(ik_psig); const int npw = this->pw_wfc_->npwk[ik]; const int total_lm = (this->p_ucell_->lmax + 1) * (this->p_ucell_->lmax + 1); ModuleBase::matrix ylm(total_lm, npw); diff --git a/source/module_psi/psi_initializer_nao.h b/source/module_psi/psi_initializer_nao.h index c0a26c2bbf..5ef7d1b6ca 100644 --- a/source/module_psi/psi_initializer_nao.h +++ b/source/module_psi/psi_initializer_nao.h @@ -15,7 +15,7 @@ class psi_initializer_nao : public psi_initializer psi_initializer_nao() {this->set_method("nao");}; ~psi_initializer_nao() {}; - virtual void proj_ao_onkG(int ik) override; + virtual void proj_ao_onkG(const int ik) override; #ifdef __MPI // MPI additional implementation /// @brief initialize the psi_initializer with external data and methods diff --git a/source/module_psi/psi_initializer_nao_random.cpp b/source/module_psi/psi_initializer_nao_random.cpp index 25f163d060..15f376fcd5 100644 --- a/source/module_psi/psi_initializer_nao_random.cpp +++ b/source/module_psi/psi_initializer_nao_random.cpp @@ -11,9 +11,6 @@ void psi_initializer_nao_random::initialize(Structure_Factor* sf, const int& rank) { psi_initializer_nao::initialize(sf, pw_wfc, p_ucell, p_parakpts, random_seed, p_pspot_nl, rank); - this->ixy2is_.clear(); - this->ixy2is_.resize(this->pw_wfc_->fftnxy); - this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); } #else template @@ -24,17 +21,15 @@ void psi_initializer_nao_random::initialize(Structure_Factor* sf, pseudopot_cell_vnl* p_pspot_nl) { psi_initializer_nao::initialize(sf, pw_wfc, p_ucell, random_seed, p_pspot_nl); - this->ixy2is_.clear(); - this->ixy2is_.resize(this->pw_wfc_->fftnxy); - this->pw_wfc_->getfftixy2is(this->ixy2is_.data()); } #endif template -void psi_initializer_nao_random::proj_ao_onkG(int ik) +void psi_initializer_nao_random::proj_ao_onkG(const int ik) { double rm = this->random_mix(); - this->psig_->fix_k(ik); + const int ik_psig = (this->psig_->get_nk() == 1) ? 0 : ik; + this->psig_->fix_k(ik_psig); psi_initializer_nao::proj_ao_onkG(ik); psi::Psi psi_random(1, this->psig_->get_nbands(), this->psig_->get_nbasis(), nullptr); psi_random.fix_k(0); diff --git a/source/module_psi/psi_initializer_nao_random.h b/source/module_psi/psi_initializer_nao_random.h index 736d49c5d2..8fd6edb4c6 100644 --- a/source/module_psi/psi_initializer_nao_random.h +++ b/source/module_psi/psi_initializer_nao_random.h @@ -34,7 +34,7 @@ class psi_initializer_nao_random : public psi_initializer_nao pseudopot_cell_vnl* = nullptr) override;//< nonlocal pseudopotential #endif - virtual void proj_ao_onkG(int ik) override; + virtual void proj_ao_onkG(const int ik) override; virtual void tabulate() override {psi_initializer_nao::tabulate();}; }; #endif \ No newline at end of file diff --git a/source/module_psi/psi_initializer_random.cpp b/source/module_psi/psi_initializer_random.cpp index c02d3b7607..fcabca18fc 100644 --- a/source/module_psi/psi_initializer_random.cpp +++ b/source/module_psi/psi_initializer_random.cpp @@ -56,10 +56,11 @@ void psi_initializer_random::random(T* psi, } template -void psi_initializer_random::proj_ao_onkG(int ik) +void psi_initializer_random::proj_ao_onkG(const int ik) { ModuleBase::timer::tick("psi_initializer_random", "initialize"); - this->psig_->fix_k(ik); + const int ik_psig = (this->psig_->get_nk() == 1) ? 0 : ik; + this->psig_->fix_k(ik_psig); this->random(this->psig_->get_pointer(), 0, this->psig_->get_nbands(), ik); ModuleBase::timer::tick("psi_initializer_random", "initialize"); } diff --git a/source/module_psi/psi_initializer_random.h b/source/module_psi/psi_initializer_random.h index 0c1f53edbd..39707a303c 100644 --- a/source/module_psi/psi_initializer_random.h +++ b/source/module_psi/psi_initializer_random.h @@ -23,7 +23,7 @@ class psi_initializer_random : public psi_initializer /// @brief calculate and output planewave wavefunction /// @param ik kpoint index /// @return initialized planewave wavefunction (psi::Psi>*) - virtual void proj_ao_onkG(int ik) override; + virtual void proj_ao_onkG(const int ik) override; #ifdef __MPI // MPI additional implementation /// @brief initialize the psi_initializer with external data and methods virtual void initialize(Structure_Factor*, //< structure factor From 20d283b14d7949a027981cc9cfd7ec8ac3ce0f29 Mon Sep 17 00:00:00 2001 From: liiutao <74701833+A-006@users.noreply.github.com> Date: Thu, 19 Sep 2024 11:10:31 +0800 Subject: [PATCH 7/8] move back Kpoint when basis_type="pw" (#5129) * move back gamma_only and Kpoint * add gtest in the OFS --- source/module_io/read_input_item_elec_stru.cpp | 17 ++++++++++++++--- source/module_io/test_serial/CMakeLists.txt | 1 - .../test_serial/read_input_item_test.cpp | 7 +++++++ 3 files changed, 21 insertions(+), 4 deletions(-) diff --git a/source/module_io/read_input_item_elec_stru.cpp b/source/module_io/read_input_item_elec_stru.cpp index cd6f87e2c4..4e7795a853 100644 --- a/source/module_io/read_input_item_elec_stru.cpp +++ b/source/module_io/read_input_item_elec_stru.cpp @@ -451,10 +451,21 @@ void ReadInput::item_elec_stru() "set to 1, a fast algorithm is used"; read_sync_bool(input.gamma_only); item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.gamma_only && para.input.basis_type == "pw") + if (para.input.basis_type == "pw" && para.input.gamma_only) { - para.input.gamma_only = false; - GlobalV::ofs_warning << "gamma_only is not supported in the pw model" << std::endl; + para.input.gamma_only = false; + GlobalV::ofs_warning << " WARNING : gamma_only has not been implemented for pw yet" << std::endl; + GlobalV::ofs_warning << "gamma_only is not supported in the pw model" << std::endl; + GlobalV::ofs_warning << " the INPUT parameter gamma_only has been reset to 0" << std::endl; + GlobalV::ofs_warning << " and a new KPT is generated with " + "gamma point as the only k point"<< std::endl; + GlobalV::ofs_warning << " Auto generating k-points file: " << para.input.kpoint_file << std::endl; + std::ofstream ofs(para.input.kpoint_file.c_str()); + ofs << "K_POINTS" << std::endl; + ofs << "0" << std::endl; + ofs << "Gamma" << std::endl; + ofs << "1 1 1 0 0 0" << std::endl; + ofs.close(); } }; this->add_item(item); diff --git a/source/module_io/test_serial/CMakeLists.txt b/source/module_io/test_serial/CMakeLists.txt index ac0a04b1b9..e6a7f3b80d 100644 --- a/source/module_io/test_serial/CMakeLists.txt +++ b/source/module_io/test_serial/CMakeLists.txt @@ -33,7 +33,6 @@ AddTest( AddTest( TARGET io_read_item_serial LIBS parameter ${math_libs} base device io_input_serial - LIBS parameter ${math_libs} base device io_input_serial SOURCES read_input_item_test.cpp ) diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp index e0eafcf0d3..9f7c6abe89 100644 --- a/source/module_io/test_serial/read_input_item_test.cpp +++ b/source/module_io/test_serial/read_input_item_test.cpp @@ -811,8 +811,15 @@ TEST_F(InputTest, Item_test) auto it = find_label("gamma_only", readinput.input_lists); param.input.basis_type = "pw"; param.input.gamma_only = true; + std::string filename = param.input.kpoint_file; it->second.reset_value(it->second, param); EXPECT_EQ(param.input.gamma_only, false); + std::ifstream ifs(filename); + std::string line; + std::getline(ifs, line); + ifs.close(); + EXPECT_EQ(line, "K_POINTS"); + } { // out_mat_r auto it = find_label("out_mat_r", readinput.input_lists); From c7349b89bbc100380c3fecf52f67455c264b2050 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Thu, 19 Sep 2024 12:03:25 +0800 Subject: [PATCH 8/8] Delley's grid for quadrature on the unit sphere (#5131) * init commit * first version * fix a few typos and add a test * finish unit test * fix compilation issue * minor clean up --- source/module_base/CMakeLists.txt | 1 + source/module_base/grid/delley.cpp | 376 +++++++++++++++++++ source/module_base/grid/delley.h | 57 +++ source/module_base/grid/test/CMakeLists.txt | 9 + source/module_base/grid/test/test_delley.cpp | 138 +++++++ source/module_base/ylm.cpp | 6 +- 6 files changed, 584 insertions(+), 3 deletions(-) create mode 100644 source/module_base/grid/delley.cpp create mode 100644 source/module_base/grid/delley.h create mode 100644 source/module_base/grid/test/CMakeLists.txt create mode 100644 source/module_base/grid/test/test_delley.cpp diff --git a/source/module_base/CMakeLists.txt b/source/module_base/CMakeLists.txt index 065c2cbe3c..dde98f915f 100644 --- a/source/module_base/CMakeLists.txt +++ b/source/module_base/CMakeLists.txt @@ -73,6 +73,7 @@ if(BUILD_TESTING) add_subdirectory(kernels/test) add_subdirectory(module_mixing/test) add_subdirectory(module_device/test) + add_subdirectory(grid/test) if (USE_ABACUS_LIBM) add_subdirectory(libm/test) endif() diff --git a/source/module_base/grid/delley.cpp b/source/module_base/grid/delley.cpp new file mode 100644 index 0000000000..1e0b6c703e --- /dev/null +++ b/source/module_base/grid/delley.cpp @@ -0,0 +1,376 @@ +#include "module_base/grid/delley.h" + +#include +#include +#include + +namespace { + +struct Table { + const int lmax_; + const int ngrid_; + const int ntype_[6]; + const std::vector data_; +}; + +// Delley's table from the original article +const std::vector table = { + { + 17, 110, {1, 1, 0, 3, 1, 0}, + { + 0.00000000000000000, 0.00000000000000000, 0.0038282704949371616, + 0.57735026918962576, 0.57735026918962576, 0.0097937375124875125, + 0.18511563534473617, 0.18511563534473617, 0.0082117372831911110, + 0.39568947305594191, 0.39568947305594191, 0.0095954713360709628, + 0.69042104838229218, 0.21595729184584883, 0.0099428148911781033, + 0.47836902881215020, 0.00000000000000000, 0.0096949963616630283, + } + }, + { + 23, 194, {1, 1, 1, 4, 1, 1}, + { + 0.00000000000000000, 0.00000000000000000, 0.0017823404472446112, + 0.57735026918962576, 0.57735026918962576, 0.0055733831788487380, + 0.70710678118654752, 0.00000000000000000, 0.0057169059499771019, + 0.44469331787174373, 0.44469331787174373, 0.0055187714672736137, + 0.28924656275754386, 0.28924656275754386, 0.0051582377118053831, + 0.67129734426952263, 0.31419699418258608, 0.0056087040825879968, + 0.12993354476500669, 0.12993354476500669, 0.0041067770281693941, + 0.34577021976112827, 0.00000000000000000, 0.0050518460646148085, + 0.52511857244364202, 0.15904171053835295, 0.0055302489162330937, + } + }, + { + 29, 302, {1, 1, 0, 6, 2, 2}, + { + 0.00000000000000000, 0.00000000000000000, 0.0008545911725128148, + 0.57735026918962576, 0.57735026918962576, 0.0035991192850255715, + 0.70117664160895449, 0.12923867271051493, 0.0036500458076772554, + 0.65663294102196118, 0.37103417838482119, 0.0036048226014198817, + 0.47290541325810046, 0.47290541325810046, 0.0035767296617433671, + 0.35156403455701051, 0.35156403455701051, 0.0034497884243058833, + 0.22196452362941784, 0.22196452362941784, 0.0031089531224136753, + 0.09618308522614784, 0.09618308522614784, 0.0023521014136891644, + 0.57189558918789607, 0.00000000000000000, 0.0036008209322164603, + 0.26441528870606625, 0.00000000000000000, 0.0029823449631718039, + 0.54486773725807738, 0.25100347517704651, 0.0035715405542733871, + 0.41277240831685310, 0.12335485325833274, 0.0033923122050061702, + } + }, + { + 35, 434, {1, 1, 1, 7, 2, 4}, + { + 0.00000000000000000, 0.00000000000000000, 0.0005265897968224436, + 0.57735026918962576, 0.57735026918962576, 0.0025123174189273072, + 0.70710678118654752, 0.00000000000600000, 0.0025482199720026072, + 0.69093463075091106, 0.21264682470755207, 0.0025304038011863550, + 0.64566647074242561, 0.40771266489776951, 0.0025132671745975644, + 0.49143426377847465, 0.49143426377847465, 0.0025017251684029361, + 0.39272597633680022, 0.39272597633680022, 0.0024453734373129800, + 0.28612890103076384, 0.28612890103076384, 0.0023026947822274158, + 0.17748360546091578, 0.17748360546091578, 0.0020142790209185282, + 0.07568084367178018, 0.07568084367178018, 0.0014624956215946138, + 0.21027252285730696, 0.00000000000000000, 0.0019109512821795323, + 0.47159869115131592, 0.00000000000000000, 0.0024174423756389808, + 0.33443631453434549, 0.09921769636429237, 0.0022366077604378487, + 0.45023303825826254, 0.20548236964030437, 0.0024169300443247753, + 0.55501523610768072, 0.31042840351665415, 0.0024966440545530860, + 0.59051570489252711, 0.10680182607580483, 0.0025122368545634951, + } + }, + { + 41, 590, {1, 1, 0, 8, 4, 6}, + { + 0.00000000000000000, 0.00000000000000000, 0.0001009005753378758, + 0.57735026918962576, 0.57735026918962576, 0.0018514016873890461, + 0.70404760433146996, 0.09291900596883211, 0.0018686219518306975, + 0.68084561988024238, 0.26999719217017240, 0.0018648696345606001, + 0.63723669159418917, 0.43342680786054810, 0.0018497643975168892, + 0.50447558060926046, 0.50447558060926046, 0.0018450277740822388, + 0.42175447334398773, 0.42175447334398773, 0.0018164174988262214, + 0.33201962086729379, 0.33201962086729379, 0.0017449464690023229, + 0.23917494336556047, 0.23917494336556047, 0.0016278016126848035, + 0.14024070738935403, 0.14024070738935403, 0.0015576827519901693, + 0.09161634328605240, 0.00000000000000000, 0.0012680968886048433, + 0.20326292518419433, 0.00000000000000000, 0.0011183965414769017, + 0.39364042372978295, 0.00000000000000000, 0.0017287035120530033, + 0.61262355812929648, 0.00000000000000000, 0.0018551905629473527, + 0.28114771623428322, 0.08959875911893791, 0.0014697353123693616, + 0.38175470908581117, 0.17327600238498666, 0.0016819651914742022, + 0.47452376478986998, 0.26422260656245780, 0.0017876372876796954, + 0.56127905075920534, 0.35189965873835832, 0.0018400735685528423, + 0.50324791996964975, 0.08886791018186295, 0.0018072536817113700, + 0.59768324320748616, 0.18154345643517542, 0.0018527289739424312, + } + }, + { + 47, 770, {1, 1, 1, 9, 4, 9}, + { + 0.00000000000000000, 0.00000000000000000, 0.0011685335608691628, + 0.57735026918962576, 0.57735026918962576, 0.0014121215930643264, + 0.70710678118654752, 0.00000000000000000, 0.0014468645950992776, + 0.11441365123336336, 0.11441365123336336, 0.0010478418864629224, + 0.19944675708548970, 0.19944675708548970, 0.0012392547584848484, + 0.28401278368259530, 0.28401278368259530, 0.0013259295792415379, + 0.36646411416548296, 0.36646411416548296, 0.0013756097758625958, + 0.44356118052513995, 0.44356118052513995, 0.0013999348863558624, + 0.51435709575333968, 0.51435709575333968, 0.0014096221218822673, + 0.63052081196671812, 0.45264446462279973, 0.0014108746499638577, + 0.67164784337293865, 0.31269529735024947, 0.0014134887639034478, + 0.69812332010174177, 0.15889512220405632, 0.0014366946685816802, + 0.12047667931264991, 0.00000000000000000, 0.0010901543574180667, + 0.30940302315480606, 0.00000000000000000, 0.0001869137844803852, + 0.34884276430183016, 0.00000000000000000, 0.0011284267652336505, + 0.53224214285417946, 0.00000000000000000, 0.0013844558026568455, + 0.23249923409267532, 0.06616159933437003, 0.0011853923885095502, + 0.32477344409682044, 0.14568618765136356, 0.0012949021664637693, + 0.41056989039349425, 0.22832839132127622, 0.0013525857420363760, + 0.49213658085114203, 0.30714431901543855, 0.0013925025908786082, + 0.56548849812588755, 0.38271180625074657, 0.0014073257894372725, + 0.43713473693946563, 0.07970715187939190, 0.0013128954307755017, + 0.52320749473197761, 0.15892620239864833, 0.0013784632898490457, + 0.60283033994386521, 0.23667220253873893, 0.0014125450609821936, + 0.62037164721742807, 0.07982328826030880, 0.0014289835314095131, + } + }, + { + 53, 974, {1, 1, 0, 12, 4, 12}, + { + 0.00000000000000000, 0.00000000800000000, 0.0001438294190527431, + 0.57735026918962576, 0.57735026918962576, 0.0011257722882870041, + 0.04292963545341347, 0.04292963545341347, 0.0004948029341949241, + 0.10514268540864042, 0.10514268540864042, 0.0007357990109125470, + 0.17500248676230874, 0.17500248676230874, 0.0008889132771304384, + 0.24776533796502568, 0.24776533796502568, 0.0009888347838921435, + 0.32065671239559574, 0.32065671239559574, 0.0010532996817094706, + 0.39165207498499835, 0.39165207498499835, 0.0010927788070145785, + 0.45908258741876237, 0.45908258741876237, 0.0011143893940632272, + 0.52145638884158605, 0.52145638884158605, 0.0011237247880515553, + 0.62531702446541989, 0.46685890569574328, 0.0011252393252438136, + 0.66379267445231699, 0.34461365423743822, 0.0011261532718159050, + 0.69104103984983007, 0.21195415185018465, 0.0011302869311238408, + 0.70529070074577603, 0.07162440144995566, 0.0011349865343639549, + 0.12366867626579899, 0.00000000000000000, 0.0006823367927109931, + 0.29407771144683870, 0.00000000000000000, 0.0009454158160447096, + 0.46977538492076491, 0.00000000000000000, 0.0010744299753856791, + 0.63345632411395669, 0.00000000000000000, 0.0011293000865691317, + 0.20291287527775228, 0.05974048614181342, 0.0008436884500901954, + 0.46026219424840539, 0.13757604084736365, 0.0010752557204488846, + 0.50306739996620357, 0.33910165263362857, 0.0011085772368644620, + 0.28176064224421343, 0.12716751914398195, 0.0009566475323783357, + 0.43315612917201574, 0.26931207404135125, 0.0010806632507173907, + 0.62561673585808142, 0.14197864526019183, 0.0011267971311962946, + 0.37983952168591567, 0.06709284600738255, 0.0010225687153580612, + 0.55175054214235205, 0.07057738183256172, 0.0011089602677131075, + 0.60296191561591869, 0.27838884778821546, 0.0011227906534357658, + 0.35896063295890958, 0.19795789389174069, 0.0010324018471174598, + 0.53486664381354765, 0.20873070611032740, 0.0011072493822838539, + 0.56749975460743735, 0.40551221378728359, 0.0011217800485199721, + } + }, + { + 59, 1202, {1, 1, 1, 13, 4, 16}, + { + 0.00006000000000000, 0.00000000000000000, 0.0001105189233267572, + 0.57735026918962576, 0.57735026918962576, 0.0009133159786443561, + 0.70710678118654752, 0.00000000000000000, 0.0009205232738090741, + 0.03712636449657089, 0.03712636449657089, 0.0003690421898017899, + 0.09140060412262223, 0.09140060412262223, 0.0005603990928680660, + 0.15310778524699062, 0.15310778524699062, 0.0006865297629282609, + 0.21809288916606116, 0.21809288916606116, 0.0007720338551145630, + 0.28398745322001746, 0.28398745322001746, 0.0008301545958894795, + 0.34911776009637644, 0.34911776009637644, 0.0008686692550179628, + 0.41214314614443092, 0.41214314614443092, 0.0008927076285846890, + 0.47189936271491266, 0.47189936271491266, 0.0009060820238568219, + 0.52731454528423366, 0.52731454528423366, 0.0009119777254940867, + 0.62094753324440192, 0.47838093807695216, 0.0009128720138604181, + 0.65697227118572905, 0.36983086645942597, 0.0009130714935691735, + 0.68417883090701434, 0.25258395570071777, 0.0009152873784554116, + 0.70126043301236308, 0.12832618665972300, 0.0009187436274321654, + 0.10723822154781661, 0.00000000000000000, 0.0005176977312965694, + 0.25820689594969680, 0.00006000000000000, 0.0007331143682101417, + 0.41727529553067168, 0.00000000000000000, 0.0008463232836379928, + 0.57003669117925033, 0.00000000000000000, 0.0009031122694253992, + 0.17717740226153253, 0.05210639477011284, 0.0006485778453163257, + 0.24757164634262876, 0.11156409571564867, 0.0007435030910982369, + 0.31736152466119767, 0.17465516775786261, 0.0008101731497468018, + 0.38542911506692237, 0.23902784793817240, 0.0008556299257311812, + 0.45074225931570644, 0.30294669735289819, 0.0008850282341265444, + 0.51235184864198708, 0.36498322605976536, 0.0009022692938426915, + 0.56937024984684411, 0.42386447815223403, 0.0009105760258970126, + 0.33546162890664885, 0.05905888853235508, 0.0007998527891839054, + 0.40902684270853572, 0.12172350510959870, 0.0008483389574594331, + 0.47853206759224352, 0.18575051945473351, 0.0008811048182425720, + 0.54343035696939004, 0.24941121623622365, 0.0009010091677105086, + 0.60311616930963100, 0.31122759471496082, 0.0009107813579482705, + 0.49322211848512846, 0.06266250624154169, 0.0008803208679738260, + 0.56321230207620997, 0.12677748006842827, 0.0009021342299040653, + 0.62698055090243917, 0.19060182227792370, 0.0009131578003189435, + 0.63942796347491023, 0.06424549224220589, 0.0009158016174693465, + } + } +}; // end of the definition of "table" + +// size of each group of points with octahedral symmetry +// 6: (1, 0, 0) x sign x permutation (vertices) +// 8: (sqrt(1/3), sqrt(1/3), sqrt(1/3)) x sign (face centers) +// 12: (sqrt(2)/2, sqrt(2)/2, 0) x sign x permutation (edge centers) +// 24a: (u, u, sqrt(1-2u^2)) x sign x permutation +// 24b: (u, 0, sqrt(1-u^2)) x sign x permutation +// 48: (u, v, sqrt(1-u^2-v^2)) x sign x permutation +const int group_size[] = {6, 8, 12, 24, 24, 48}; + +using Fill_t = std::function; + +// functors that fill the grid group-wise +const std::vector fill = { + // (1, 0, 0) x sign x permutation (vertices) + [](double* grid, double, double) { + for (int i = 0; i < 3; ++i) { + for (double one : {-1.0, 1.0}) { + grid[i] = one; + grid[(i+1)%3] = 0.0; + grid[(i+2)%3] = 0.0; + std::advance(grid, 3); + } + } + }, + + // (sqrt(1/3), sqrt(1/3), sqrt(1/3)) x sign (face centers) + [](double* grid, double, double) { + const double a = std::sqrt(3) / 3.0; + for (int xsign : {-1, 1}) { + for (int ysign : {-1, 1}) { + for (int zsign : {-1, 1}) { + grid[0] = xsign * a; + grid[1] = ysign * a; + grid[2] = zsign * a; + std::advance(grid, 3); + } + } + } + }, + + // (sqrt(2)/2, sqrt(2)/2, 0) x sign x permutation (edge centers) + [](double* grid, double, double) { + const double a = std::sqrt(2) / 2.0; + for (int i = 0; i < 3; ++i) { + for (int sign1 : {-1, 1}) { + for (int sign2 : {-1, 1}) { + grid[i] = 0; + grid[(i+1)%3] = sign1 * a; + grid[(i+2)%3] = sign2 * a; + std::advance(grid, 3); + } + } + } + }, + + // (u, u, sqrt(1-2u^2)) x sign x permutation + [](double* grid, double x, double y) { + double u = x == y ? x : std::sqrt(1.0 - x * x - y * y); + double v = std::sqrt(1.0 - 2.0 * u * u); + for (int i = 0; i < 3; ++i) { + for (int sign1 : {-1, 1}) { + for (int sign2 : {-1, 1}) { + for (int sign3 : {-1, 1}) { + grid[i] = sign1 * u; + grid[(i+1)%3] = sign2 * u; + grid[(i+2)%3] = sign3 * v; + std::advance(grid, 3); + } + } + } + } + }, + + // (u, 0, sqrt(1-u^2)) x sign x permutation + [](double* grid, double x, double y) { + double u = x > 0 ? x : y; + double v = std::sqrt(1.0 - u * u); + for (int i0 = 0; i0 < 3; ++i0) { + for (int iu0 : {1, 2}) { + for (int sign_u : {-1, 1}) { + for (int sign_v : {-1, 1}) { + grid[i0] = 0; + grid[(i0+iu0)%3] = sign_u * u; + grid[(i0-iu0+3)%3] = sign_v * v; + std::advance(grid, 3); + } + } + } + } + }, + + // (u, v, sqrt(1-u^2-v^2)) x sign x permutation + [](double* grid, double x, double y) { + double r = x; + double s = y; + double t = std::sqrt(1.0 - r * r - s * s); + for (int ir = 0; ir < 3; ++ir) { + for (int irs : {1, 2}) { + for (int sign_r : {-1, 1}) { + for (int sign_s : {-1, 1}) { + for (int sign_t : {-1, 1}) { + grid[ir] = sign_r * r; + grid[(ir+irs)%3] = sign_s * s; + grid[(ir-irs+3)%3] = sign_t * t; + std::advance(grid, 3); + } + } + } + } + } + }, +}; // end of the definition of "fill" + +const Table* _find_table(int& lmax) { + // NOTE: this function assumes elements in "Delley::table_" are + // arranged such that their members "lmax_" are in ascending order. + auto tab = std::find_if(table.begin(), table.end(), + [lmax](const Table& t) { return t.lmax_ >= lmax; }); + return tab == table.end() ? nullptr : (lmax = tab->lmax_, &(*tab)); +} + +void _get(const Table* tab, double* grid, double* weight) { + assert(tab); + const double* ptr = &tab->data_[0]; + for (int itype = 0; itype < 6; ++itype) { + int stride = group_size[itype]; + for (int i = 0; i < tab->ntype_[itype]; ++i, ptr += 3, + grid += 3*stride, weight += stride) { + fill[itype](grid, ptr[0], ptr[1]); + std::fill(weight, weight + stride, ptr[2]); + } + } +} + +} // end of anonymous namespace + + +int Grid::Delley::ngrid(int& lmax) { + auto tab = _find_table(lmax); + return tab ? tab->ngrid_ : -1; +} + + +int Grid::Delley::get(int& lmax, double* grid, double* weight) { + auto tab = _find_table(lmax); + return tab ? _get(tab, grid, weight), 0 : -1; +} + + +int Grid::Delley::get( + int& lmax, + std::vector& grid, + std::vector& weight +) { + auto tab = _find_table(lmax); + if (!tab) { + return -1; + } + grid.resize(3 * tab->ngrid_); + weight.resize(tab->ngrid_); + _get(tab, grid.data(), weight.data()); + return 0; +} diff --git a/source/module_base/grid/delley.h b/source/module_base/grid/delley.h new file mode 100644 index 0000000000..4f68a2168a --- /dev/null +++ b/source/module_base/grid/delley.h @@ -0,0 +1,57 @@ +#ifndef GRID_DELLEY_H +#define GRID_DELLEY_H + +#include +#include + +/** + * @brief Delley's grid for quadrature on the unit sphere. + * + * Reference: + * Delley, B. (1996). High order integration schemes on the unit sphere. + * Journal of computational chemistry, 17(9), 1152-1155. + * + */ +namespace Grid { +namespace Delley { + +/** + * @brief Number of Delley's grid points for a certain order of accuracy. + * + * This function finds the minimum Delley's grid with an accuracy + * of order "lmax". On exit, lmax is set to the order of this grid, and + * its corresponding number of grid points is returned. If no such grid + * is available, lmax is left unchanged and the function will return -1. + * + * For example, if lmax = 20 on input, the function will return 194 and + * lmax will be set to 23. + * + */ +int ngrid(int& lmax); + + +/** + * @brief Delley's quadrature grid and weights. + * + * This function retrieves the minimum Delley's grid with an accuray + * of order "lmax". On exit, lmax is set to the order of this grid, and + * the coordinates & weights are returned in "grid" & "weight". + * + * Coordinates are stored in the following order: + * + * x0, y0, z0, x1, y1, z1, x2, y2, z2, ... + * + * "grid" and "weight" must be pre-allocated to hold 3*ngrid(lmax) and + * ngrid(lmax) elements, respectively. The function will return 0 if + * successful, or -1 if the requested order cannot be fulfilled. + * + */ +int get(int& lmax, double* grid, double* weight); + +// a handy wrapper doing the same as above +int get(int& lmax, std::vector& grid, std::vector& weight); + +} // end of namespace Delley +} // end of namespace Grid + +#endif diff --git a/source/module_base/grid/test/CMakeLists.txt b/source/module_base/grid/test/CMakeLists.txt new file mode 100644 index 0000000000..0f77a425fe --- /dev/null +++ b/source/module_base/grid/test/CMakeLists.txt @@ -0,0 +1,9 @@ +remove_definitions(-D__MPI) + +AddTest( + TARGET test_delley + SOURCES test_delley.cpp + ../delley.cpp + ../../ylm.cpp +) + diff --git a/source/module_base/grid/test/test_delley.cpp b/source/module_base/grid/test/test_delley.cpp new file mode 100644 index 0000000000..2b892cd532 --- /dev/null +++ b/source/module_base/grid/test/test_delley.cpp @@ -0,0 +1,138 @@ +#include "module_base/grid/delley.h" +#include "module_base/ylm.h" + +#include "gtest/gtest.h" +#include +#ifdef __MPI +#include +#endif + +using namespace Grid; + +// mock the function to prevent unnecessary dependency +namespace ModuleBase { +void WARNING_QUIT(const std::string&, const std::string&) {} +} + +class DelleyTest: public ::testing::Test { +protected: + void randgen(int lmax, std::vector& coef); + const double tol = 1e-12; +}; + + +void DelleyTest::randgen(int lmax, std::vector& coef) { + coef.resize((lmax + 1) * (lmax + 1)); + + // fill coef with uniformly distributed random numbers + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution dis(0.0, 1.0); + for (size_t i = 0; i < coef.size(); ++i) { + coef[i] = dis(gen); + } + + // normalize the coefficients + double fac = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + fac += coef[i] * coef[i]; + } + + fac = 1.0 / std::sqrt(fac); + for (size_t i = 0; i < coef.size(); ++i) { + coef[i] *= fac; + } +} + + +TEST_F(DelleyTest, NumGrid) { + int lmax = 5; + int ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 17); + EXPECT_EQ(ngrid, 110); + + lmax = 17; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 17); + EXPECT_EQ(ngrid, 110); + + lmax = 20; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 23); + EXPECT_EQ(ngrid, 194); + + lmax = 59; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 59); + EXPECT_EQ(ngrid, 1202); + + lmax = 60; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 60); + EXPECT_EQ(ngrid, -1); +} + + +TEST_F(DelleyTest, Accuracy) { + /* + * Given + * + * f = c[0]*Y00 + c[1]*Y10 + c[2]*Y11 + ..., + * + * where c[0], c[1], c[2], ... are some random numbers, the integration + * of |f|^2 on the unit sphere + * + * \int |f|^2 d\Omega = c[0]^2 + c[1]^2 + c[2]^2 + ... . + * + * This test verifies with the above integral that quadrature with + * Delley's grid is exact up to floating point errors. + * + */ + std::vector grid, weight, coef; + + for (int grid_lmax = 17; grid_lmax < 60; grid_lmax +=6) { + Delley::get(grid_lmax, grid, weight); + int func_lmax = grid_lmax / 2; + randgen(func_lmax, coef); + + double val = 0.0; + std::vector ylm_real; + for (size_t i = 0; i < weight.size(); i++) { + ModuleBase::Ylm::sph_harm(func_lmax, + grid[3*i], grid[3*i+1], grid[3*i+2], ylm_real); + double tmp = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + tmp += coef[i] * ylm_real[i]; + } + val += weight[i] * tmp * tmp; + } + val *= 4.0 * std::acos(-1.0); + + double val_ref = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + val_ref += coef[i] * coef[i]; + } + + double abs_diff = std::abs(val - val_ref); + EXPECT_LT(abs_diff, tol); + //printf("order = %2i val_ref = %8.5f abs_diff = %8.5e\n", + // grid_lmax, val_ref, abs_diff); + } +} + + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + + return result; +} diff --git a/source/module_base/ylm.cpp b/source/module_base/ylm.cpp index a60ea9b11c..a36b14c60a 100644 --- a/source/module_base/ylm.cpp +++ b/source/module_base/ylm.cpp @@ -54,7 +54,7 @@ std::vector Ylm::ylmcoef = { // here Lmax == max angular momentum + 1 void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3 &vec, double ylmr[] ) { - ModuleBase::timer::tick ("Ylm","get_ylm_real"); + //ModuleBase::timer::tick ("Ylm","get_ylm_real"); //1e-9 is too large const double cut0 = 1e-12; // allocate space. @@ -165,7 +165,7 @@ void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3 &vec, } }// end do - ModuleBase::timer::tick ("Ylm", "get_ylm_real"); + //ModuleBase::timer::tick ("Ylm", "get_ylm_real"); return; } @@ -536,7 +536,7 @@ void Ylm::sph_harm std::vector &rly ) { - rly.reserve( (Lmax+1)*(Lmax+1) ); + rly.resize( (Lmax+1)*(Lmax+1) ); //begin calculation /***************************