diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index f6f90cb382..33be890791 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -219,7 +219,8 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p this->kv.ngk.data(), this->pw_wfc->npwk_max, &this->sf, - &this->ppcell); + &this->ppcell, + ucell); this->kspw_psi = PARAM.inp.device == "gpu" || PARAM.inp.precision == "single" ? new psi::Psi(this->psi[0]) @@ -257,7 +258,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->pw_wfc->collect_local_pw(PARAM.inp.erf_ecut, PARAM.inp.erf_height, PARAM.inp.erf_sigma); - this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell); + this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell,ucell); } if (ucell.ionic_position_updated) { @@ -373,6 +374,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->kspw_psi, this->p_hamilt, this->ppcell, + ucell, GlobalV::ofs_running, this->already_initpsi); diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp index 975525717e..770a4184e2 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.cpp @@ -61,7 +61,7 @@ void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const Modul ModuleBase::TITLE("PW_Basis","setup_structure_factor"); ModuleBase::timer::tick("PW_Basis","setup_struc_factor"); const std::complex ci_tpi = ModuleBase::NEG_IMAG_UNIT * ModuleBase::TWO_PI; - + this->ucell = Ucell; this->strucFac.create(Ucell->ntype, rho_basis->npw); ModuleBase::Memory::record("SF::strucFac", sizeof(std::complex) * Ucell->ntype*rho_basis->npw); diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h index 46ebbc1608..7018e320b7 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h +++ b/source/module_hamilt_pw/hamilt_pwdft/structure_factor.h @@ -53,6 +53,7 @@ class Structure_Factor ModuleBase::Vector3 q); private: + const UnitCell* ucell; std::complex * c_eigts1 = nullptr, * c_eigts2 = nullptr, * c_eigts3 = nullptr; std::complex * z_eigts1 = nullptr, * z_eigts2 = nullptr, * z_eigts3 = nullptr; const ModulePW::PW_Basis* rho_basis = nullptr; diff --git a/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp b/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp index 83b3f37ce0..add76f6fb3 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/structure_factor_k.cpp @@ -10,7 +10,7 @@ std::complex* Structure_Factor::get_sk(const int ik, const ModulePW::PW_Basis_K* wfc_basis) const { ModuleBase::timer::tick("Structure_Factor", "get_sk"); - const double arg = (wfc_basis->kvec_c[ik] * GlobalC::ucell.atoms[it].tau[ia]) * ModuleBase::TWO_PI; + const double arg = (wfc_basis->kvec_c[ik] * ucell->atoms[it].tau[ia]) * ModuleBase::TWO_PI; const std::complex kphase = std::complex(cos(arg), -sin(arg)); const int npw = wfc_basis->npwk[ik]; std::complex *sk = new std::complex[npw]; @@ -26,19 +26,22 @@ std::complex* Structure_Factor::get_sk(const int ik, const int ixy = wfc_basis->is2fftixy[is]; int ix = ixy / wfc_basis->fftny; int iy = ixy % wfc_basis->fftny; - if (ix >= int(nx / 2) + 1) { + if (ix >= int(nx / 2) + 1) + { ix -= nx; -} - if (iy >= int(ny / 2) + 1) { + } + if (iy >= int(ny / 2) + 1) + { iy -= ny; -} - if (iz >= int(nz / 2) + 1) { + } + if (iz >= int(nz / 2) + 1) + { iz -= nz; -} + } ix += this->rho_basis->nx; iy += this->rho_basis->ny; iz += this->rho_basis->nz; - const int iat = GlobalC::ucell.itia2iat(it, ia); + const int iat = ucell->itia2iat(it, ia); sk[igl] = kphase * this->eigts1(iat, ix) * this->eigts2(iat, iy) * this->eigts3(iat, iz); } ModuleBase::timer::tick("Structure_Factor", "get_sk"); @@ -66,33 +69,33 @@ void Structure_Factor::get_sk(Device* ctx, int iat = 0, _npw = wfc_basis->npwk[ik], eigts1_nc = this->eigts1.nc, eigts2_nc = this->eigts2.nc, eigts3_nc = this->eigts3.nc; - int *igl2isz = nullptr, *is2fftixy = nullptr, *atom_na = nullptr, *h_atom_na = new int[GlobalC::ucell.ntype]; - FPTYPE *atom_tau = nullptr, *h_atom_tau = new FPTYPE[GlobalC::ucell.nat * 3], *kvec = wfc_basis->get_kvec_c_data(); + int *igl2isz = nullptr, *is2fftixy = nullptr, *atom_na = nullptr, *h_atom_na = new int[ucell->ntype]; + FPTYPE *atom_tau = nullptr, *h_atom_tau = new FPTYPE[ucell->nat * 3], *kvec = wfc_basis->get_kvec_c_data(); std::complex *eigts1 = this->get_eigts1_data(), *eigts2 = this->get_eigts2_data(), *eigts3 = this->get_eigts3_data(); - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell->ntype; it++) { - h_atom_na[it] = GlobalC::ucell.atoms[it].na; + h_atom_na[it] = ucell->atoms[it].na; } #ifdef _OPENMP #pragma omp parallel for #endif - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int iat = 0; iat < ucell->nat; iat++) { - int it = GlobalC::ucell.iat2it[iat]; - int ia = GlobalC::ucell.iat2ia[iat]; - auto *tau = reinterpret_cast(GlobalC::ucell.atoms[it].tau.data()); + int it = ucell->iat2it[iat]; + int ia = ucell->iat2ia[iat]; + auto *tau = reinterpret_cast(ucell->atoms[it].tau.data()); h_atom_tau[iat * 3 + 0] = static_cast(tau[ia * 3 + 0]); h_atom_tau[iat * 3 + 1] = static_cast(tau[ia * 3 + 1]); h_atom_tau[iat * 3 + 2] = static_cast(tau[ia * 3 + 2]); } if (device == base_device::GpuDevice) { - resmem_int_op()(ctx, atom_na, GlobalC::ucell.ntype); - syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, GlobalC::ucell.ntype); + resmem_int_op()(ctx, atom_na, ucell->ntype); + syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, ucell->ntype); - resmem_var_op()(ctx, atom_tau, GlobalC::ucell.nat * 3); - syncmem_var_op()(ctx, cpu_ctx, atom_tau, h_atom_tau, GlobalC::ucell.nat * 3); + resmem_var_op()(ctx, atom_tau, ucell->nat * 3); + syncmem_var_op()(ctx, cpu_ctx, atom_tau, h_atom_tau, ucell->nat * 3); igl2isz = wfc_basis->d_igl2isz_k; is2fftixy = wfc_basis->d_is2fftixy; @@ -107,7 +110,7 @@ void Structure_Factor::get_sk(Device* ctx, cal_sk_op()(ctx, ik, - GlobalC::ucell.ntype, + ucell->ntype, wfc_basis->nx, wfc_basis->ny, wfc_basis->nz, @@ -152,7 +155,7 @@ std::complex* Structure_Factor::get_skq(int ik, for (int ig = 0; ig < npw; ig++) { ModuleBase::Vector3 qkq = wfc_basis->getgpluskcar(ik, ig) + q; - double arg = (qkq * GlobalC::ucell.atoms[it].tau[ia]) * ModuleBase::TWO_PI; + double arg = (qkq * ucell->atoms[it].tau[ia]) * ModuleBase::TWO_PI; skq[ig] = std::complex(cos(arg), -sin(arg)); } diff --git a/source/module_hsolver/test/hsolver_pw_sup.h b/source/module_hsolver/test/hsolver_pw_sup.h index 89ee061c64..cea4b6118c 100644 --- a/source/module_hsolver/test/hsolver_pw_sup.h +++ b/source/module_hsolver/test/hsolver_pw_sup.h @@ -192,6 +192,7 @@ void diago_PAO_in_pw_k2( wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { for (int i = 0; i < wvf.size(); i++) { wvf.get_pointer()[i] = std::complex((float)i + 1, 0); @@ -207,6 +208,7 @@ void diago_PAO_in_pw_k2( wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { for (int i = 0; i < wvf.size(); i++) { wvf.get_pointer()[i] = std::complex((double)i + 1, 0); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 46697710ea..4b3b53a548 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -255,7 +255,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu // necessary steps in ESolver_KS::before_all_runners : symmetry and k-points if (ModuleSymmetry::Symmetry::symm_flag == 1) { - GlobalC::ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); + ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } this->kv.set(ucell,ucell.symm, PARAM.inp.kpoint_file, PARAM.inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running); @@ -318,12 +318,12 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu this->init_pot(chg_gs); // search adjacent atoms and init Gint - std::cout << "ucell.infoNL.get_rcutmax_Beta(): " << GlobalC::ucell.infoNL.get_rcutmax_Beta() << std::endl; + std::cout << "ucell.infoNL.get_rcutmax_Beta(): " << ucell.infoNL.get_rcutmax_Beta() << std::endl; double search_radius = -1.0; search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, PARAM.inp.out_level, orb.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), + ucell.infoNL.get_rcutmax_Beta(), PARAM.globalv.gamma_only_local); atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, @@ -341,7 +341,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu std::vector> dpsi_u; std::vector> d2psi_u; - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, orb, psi_u, dpsi_u, d2psi_u); + Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb, psi_u, dpsi_u, d2psi_u); this->gt_.set_pbc_grid(this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz, @@ -357,7 +357,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu this->pw_rho->ny, this->pw_rho->nplane, this->pw_rho->startz_current, - GlobalC::ucell, + ucell, GlobalC::GridD, dr_uniform, rcuts, diff --git a/source/module_lr/operator_casida/operator_lr_exx.cpp b/source/module_lr/operator_casida/operator_lr_exx.cpp index c3399a72b4..68ccd403b2 100644 --- a/source/module_lr/operator_casida/operator_lr_exx.cpp +++ b/source/module_lr/operator_casida/operator_lr_exx.cpp @@ -65,7 +65,7 @@ namespace LR for (auto cell : this->BvK_cells) { std::complex frac = RI::Global_Func::convert>(std::exp( - -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * GlobalC::ucell.latvec)))); + -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (this->kv.kvec_c.at(ik) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); for (int it1 = 0;it1 < ucell.ntype;++it1) for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ++ia1) for (int it2 = 0;it2 < ucell.ntype;++it2) diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index 890b03398a..52ba4fc751 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -68,7 +68,7 @@ namespace LR // 3. v_hxc = f_hxc * rho_trans ModuleBase::matrix vr_hxc(1, nrxx); //grid - this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks); + this->pot.lock()->cal_v_eff(rho_trans, ucell, vr_hxc, ispin_ks); LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1); // 4. V^{Hxc}_{\mu,\nu}=\int{dr} \phi_\mu(r) v_{Hxc}(r) \phi_\mu(r) @@ -76,7 +76,7 @@ namespace LR this->gint->get_hRGint()->set_zero(); this->gint->cal_gint(&inout_vlocal); this->hR->set_zero(); // clear hR for each bands - this->gint->transfer_pvpR(&*this->hR, &GlobalC::ucell); //grid to 2d block + this->gint->transfer_pvpR(&*this->hR, &ucell); //grid to 2d block ModuleBase::timer::tick("OperatorLRHxc", "grid_calculation"); } @@ -88,7 +88,7 @@ namespace LR elecstate::DensityMatrix, double> DM_trans_real_imag(&pmat, 1, kv.kvec_d, kv.get_nks() / nspin); DM_trans_real_imag.init_DMR(*this->hR); - hamilt::HContainer HR_real_imag(GlobalC::ucell, &this->pmat); + hamilt::HContainer HR_real_imag(ucell, &this->pmat); LR_Util::initialize_HR, double>(HR_real_imag, ucell, gd, orb_cutoff_); auto dmR_to_hR = [&, this](const char& type) -> void @@ -111,7 +111,7 @@ namespace LR // 3. v_hxc = f_hxc * rho_trans ModuleBase::matrix vr_hxc(1, nrxx); //grid - this->pot.lock()->cal_v_eff(rho_trans, GlobalC::ucell, vr_hxc, ispin_ks); + this->pot.lock()->cal_v_eff(rho_trans, ucell, vr_hxc, ispin_ks); // print_grid_nonzero(vr_hxc.c, this->poticab->nrxx, 10, "vr_hxc"); LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1); @@ -123,9 +123,9 @@ namespace LR // LR_Util::print_HR(*this->gint->get_hRGint(), this->ucell.nat, "VR(grid)"); HR_real_imag.set_zero(); - this->gint->transfer_pvpR(&HR_real_imag, &GlobalC::ucell, &GlobalC::GridD); + this->gint->transfer_pvpR(&HR_real_imag, &ucell, &GlobalC::GridD); // LR_Util::print_HR(HR_real_imag, this->ucell.nat, "VR(real, 2d)"); - LR_Util::set_HR_real_imag_part(HR_real_imag, *this->hR, GlobalC::ucell.nat, type); + LR_Util::set_HR_real_imag_part(HR_real_imag, *this->hR, ucell.nat, type); }; this->hR->set_zero(); dmR_to_hR('R'); //real diff --git a/source/module_psi/psi_init.cpp b/source/module_psi/psi_init.cpp index fdd08958ba..201539a54f 100644 --- a/source/module_psi/psi_init.cpp +++ b/source/module_psi/psi_init.cpp @@ -93,7 +93,8 @@ void PSIInit::allocate_psi(Psi>*& psi, const int* ngk, const int npwx, Structure_Factor* p_sf, - pseudopot_cell_vnl* p_ppcell) + pseudopot_cell_vnl* p_ppcell, + const UnitCell& ucell) { // allocate memory for std::complex datatype psi // New psi initializer in ABACUS, Developer's note: @@ -126,7 +127,7 @@ void PSIInit::allocate_psi(Psi>*& psi, // however, init_at_1 does not actually initialize the psi, instead, it is a // function to calculate a interpolate table saving overlap intergral or say // Spherical Bessel Transform of atomic orbitals. - this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at); + this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at); // similarly, wfcinit not really initialize any wavefunction, instead, it initialize // the mapping from ixy, the 1d flattened index of point on fft grid (x, y) plane, // to the index of "stick", composed of grid points. @@ -135,7 +136,10 @@ void PSIInit::allocate_psi(Psi>*& psi, } template -void PSIInit::make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell) +void PSIInit::make_table(const int nks, + Structure_Factor* p_sf, + pseudopot_cell_vnl* p_ppcell, + const UnitCell& ucell) { if (this->use_psiinitializer) { @@ -143,7 +147,7 @@ void PSIInit::make_table(const int nks, Structure_Factor* p_sf, pseud else // old initialization method, used in EXX calculation { this->wf_old.init_after_vc(nks); // reallocate wanf2, the planewave expansion of lcao - this->wf_old.init_at_1(p_sf, &p_ppcell->tab_at); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq + this->wf_old.init_at_1(ucell,p_sf, &p_ppcell->tab_at); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq } } @@ -152,6 +156,7 @@ void PSIInit::initialize_psi(Psi>* psi, psi::Psi* kspw_psi, hamilt::Hamilt* p_hamilt, const pseudopot_cell_vnl& nlpp, + const UnitCell& ucell, std::ofstream& ofs_running, const bool is_already_initpsi) { @@ -278,6 +283,7 @@ void PSIInit::initialize_psi(Psi>* psi, &this->wf_old, nlpp.tab_at, nlpp.lmaxkb, + ucell, p_hamilt); } } @@ -294,6 +300,7 @@ void PSIInit::initialize_psi(Psi>* psi, &this->wf_old, nlpp.tab_at, nlpp.lmaxkb, + ucell, p_hamilt); } } diff --git a/source/module_psi/psi_init.h b/source/module_psi/psi_init.h index df0ec83446..472ab4d0ff 100644 --- a/source/module_psi/psi_init.h +++ b/source/module_psi/psi_init.h @@ -36,10 +36,14 @@ class PSIInit const int* ngk, //< number of G-vectors in the current pool const int npwx, //< max number of plane waves of all pools Structure_Factor* p_sf, //< structure factor - pseudopot_cell_vnl* p_ppcell); //< nonlocal pseudopotential + pseudopot_cell_vnl* p_ppcell, //< nonlocal pseudopotential + const UnitCell& ucell); //< unit cell // make interpolate table - void make_table(const int nks, Structure_Factor* p_sf, pseudopot_cell_vnl* p_ppcell); + void make_table(const int nks, + Structure_Factor* p_sf, + pseudopot_cell_vnl* p_ppcell, + const UnitCell& ucell); //------------------------ only for psi_initializer -------------------- /** @@ -54,6 +58,7 @@ class PSIInit psi::Psi* kspw_psi, hamilt::Hamilt* p_hamilt, const pseudopot_cell_vnl& nlpp, + const UnitCell& ucell, std::ofstream& ofs_running, const bool is_already_initpsi); diff --git a/source/module_psi/wavefunc.cpp b/source/module_psi/wavefunc.cpp index 45a6e7ea18..31bdad80dd 100644 --- a/source/module_psi/wavefunc.cpp +++ b/source/module_psi/wavefunc.cpp @@ -140,7 +140,7 @@ void wavefunc::wfcinit(psi::Psi>* psi_in, ModulePW::PW_Basi return; } -int wavefunc::get_starting_nw() const +int wavefunc::get_starting_nw(const int natomwfc) const { if (PARAM.inp.init_wfc == "file") { @@ -148,7 +148,7 @@ int wavefunc::get_starting_nw() const } else if (PARAM.inp.init_wfc.substr(0, 6) == "atomic") { - if (GlobalC::ucell.natomwfc >= PARAM.inp.nbands) + if (natomwfc >= PARAM.inp.nbands) { if (PARAM.inp.test_wf) { @@ -160,11 +160,11 @@ int wavefunc::get_starting_nw() const if (PARAM.inp.test_wf) { GlobalV::ofs_running << " Start wave functions are atomic + " - << PARAM.inp.nbands - GlobalC::ucell.natomwfc << " random wave functions." + << PARAM.inp.nbands - natomwfc << " random wave functions." << std::endl; } } - return std::max(GlobalC::ucell.natomwfc, PARAM.inp.nbands); + return std::max(natomwfc, PARAM.inp.nbands); } else if (PARAM.inp.init_wfc == "random") { @@ -194,6 +194,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { // TODO float func @@ -207,6 +208,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_CPU>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -252,7 +254,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, } } else if (PARAM.inp.init_wfc == "random" - || (PARAM.inp.init_wfc.substr(0, 6) == "atomic" && GlobalC::ucell.natomwfc == 0)) + || (PARAM.inp.init_wfc.substr(0, 6) == "atomic" && ucell.natomwfc == 0)) { p_wf->random(wvf.get_pointer(), 0, nbands, ik, wfc_basis); @@ -273,7 +275,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, } else if (PARAM.inp.init_wfc.substr(0, 6) == "atomic") { - const int starting_nw = p_wf->get_starting_nw(); + const int starting_nw = p_wf->get_starting_nw(ucell.natomwfc); if (starting_nw == 0) { return; @@ -287,9 +289,9 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "starting_nw", starting_nw); } - p_wf->atomic_wfc(ik, + p_wf->atomic_wfc(ucell, + ik, current_nbasis, - GlobalC::ucell.lmax_ppwf, lmaxkb, wfc_basis, wfcatom, @@ -298,7 +300,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, PARAM.globalv.dq); if (PARAM.inp.init_wfc == "atomic+random" - && starting_nw == GlobalC::ucell.natomwfc) // added by qianrui 2021-5-16 + && starting_nw == ucell.natomwfc) // added by qianrui 2021-5-16 { p_wf->atomicrandom(wfcatom, 0, starting_nw, ik, wfc_basis); } @@ -307,7 +309,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_CPU* ctx, // If not enough atomic wfc are available, complete // with random wfcs //==================================================== - p_wf->random(wfcatom.c, GlobalC::ucell.natomwfc, nbands, ik, wfc_basis); + p_wf->random(wfcatom.c, ucell.natomwfc, nbands, ik, wfc_basis); // (7) Diago with cg method. // if(GlobalV::DIAGO_TYPE == "cg") xiaohui modify 2013-09-02 @@ -351,6 +353,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_GPU>* phm_in) { // TODO float @@ -364,6 +367,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, base_device::DEVICE_GPU>* phm_in) { ModuleBase::TITLE("wavefunc", "diago_PAO_in_pw_k2"); @@ -382,7 +386,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, ModuleIO::read_wfc_pw(filename.str(), wfc_basis, ik, p_wf->nkstot, wfcatom); } - starting_nw = p_wf->get_starting_nw(); + starting_nw = p_wf->get_starting_nw(ucell.natomwfc); if (starting_nw == 0) return; assert(starting_nw > 0); @@ -391,9 +395,9 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "starting_nw", starting_nw); if (PARAM.inp.init_wfc.substr(0, 6) == "atomic") { - p_wf->atomic_wfc(ik, + p_wf->atomic_wfc(ucell, + ik, current_nbasis, - GlobalC::ucell.lmax_ppwf, lmaxkb, wfc_basis, wfcatom, @@ -401,7 +405,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, PARAM.globalv.nqx, PARAM.globalv.dq); if (PARAM.inp.init_wfc == "atomic+random" - && starting_nw == GlobalC::ucell.natomwfc) // added by qianrui 2021-5-16 + && starting_nw == ucell.natomwfc) // added by qianrui 2021-5-16 { p_wf->atomicrandom(wfcatom, 0, starting_nw, ik, wfc_basis); } @@ -410,7 +414,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, // If not enough atomic wfc are available, complete // with random wfcs //==================================================== - p_wf->random(wfcatom.c, GlobalC::ucell.natomwfc, nbands, ik, wfc_basis); + p_wf->random(wfcatom.c, ucell.natomwfc, nbands, ik, wfc_basis); } else if (PARAM.inp.init_wfc == "random") { diff --git a/source/module_psi/wavefunc.h b/source/module_psi/wavefunc.h index adc9146bc6..2c3c21fb08 100644 --- a/source/module_psi/wavefunc.h +++ b/source/module_psi/wavefunc.h @@ -21,7 +21,7 @@ class wavefunc : public WF_atomic void wfcinit(psi::Psi>* psi_in, ModulePW::PW_Basis_K* wfc_basis); - int get_starting_nw(void) const; + int get_starting_nw(const int natomwfc) const; void init_after_vc(const int nks); // LiuXh 20180515 }; @@ -57,6 +57,7 @@ void diago_PAO_in_pw_k2(const Device* ctx, wavefunc* p_wf, const ModuleBase::realArray& tab_at, const int& lmaxkb, + const UnitCell& ucell, hamilt::Hamilt, Device>* phm_in = nullptr); } // namespace hamilt diff --git a/source/module_psi/wf_atomic.cpp b/source/module_psi/wf_atomic.cpp index fbfc9837b6..8a3b0d62a6 100644 --- a/source/module_psi/wf_atomic.cpp +++ b/source/module_psi/wf_atomic.cpp @@ -32,7 +32,9 @@ WF_atomic::~WF_atomic() } } -void WF_atomic::init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at) +void WF_atomic::init_at_1(const UnitCell& ucell, + Structure_Factor* sf_in, + ModuleBase::realArray* tab_at) { if (PARAM.inp.use_paw) { @@ -46,15 +48,15 @@ void WF_atomic::init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at this->psf = sf_in; GlobalV::ofs_running << "\n Make real space PAO into reciprocal space." << std::endl; - this->print_PAOs(); + this->print_PAOs(ucell); //---------------------------------------------------------- // EXPLAIN : Find the type of atom that has most mesh points. //---------------------------------------------------------- int ndm = 0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - ndm = (GlobalC::ucell.atoms[it].ncpp.msh > ndm) ? GlobalC::ucell.atoms[it].ncpp.msh : ndm; + ndm = (ucell.atoms[it].ncpp.msh > ndm) ? ucell.atoms[it].ncpp.msh : ndm; } ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "max mesh points in Pseudopotential", ndm); @@ -66,21 +68,21 @@ void WF_atomic::init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at // orbitals (controlled by parameters) // // USE GLOBAL CLASS VARIABLES : - // NAME : GlobalC::ucell.atoms.nchi - // NAME : GlobalC::ucell.atoms.msh(number of mesh points) - // NAME : GlobalC::ucell.atoms.r + // NAME : ucell.atoms.nchi + // NAME : ucell.atoms.msh(number of mesh points) + // NAME : ucell.atoms.r //---------------------------------------------------------- const int startq = 0; - const double pref = ModuleBase::FOUR_PI / sqrt(GlobalC::ucell.omega); + const double pref = ModuleBase::FOUR_PI / sqrt(ucell.omega); double* aux = new double[ndm]; double* vchi = new double[ndm]; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "dq(describe PAO in reciprocal space)", PARAM.globalv.dq); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "max q", PARAM.globalv.nqx); - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - Atom* atom = &GlobalC::ucell.atoms[it]; + Atom* atom = &ucell.atoms[it]; GlobalV::ofs_running << "\n number of pseudo atomic orbitals for " << atom->label << " is " << atom->ncpp.nchi << std::endl; @@ -206,13 +208,13 @@ void WF_atomic::init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at // if( it == 0 && ic == 0 ) // { // - // for (ir = 0;ir < GlobalC::ucell.atoms[it].msh;ir++) + // for (ir = 0;ir < ucell.atoms[it].msh;ir++) // GlobalV::ofs_running << std::setprecision(20) << "\n vchi(" << ir << ")=" << // vchi[ir]; GlobalV::ofs_running<<"\n aux[0] = "<* aux = new std::complex[np]; double* chiaux = nullptr; @@ -294,20 +296,20 @@ void WF_atomic::atomic_wfc(const int& ik, // Calculate G space 3D wave functions //--------------------------------------------------------- double* flq = new double[np]; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int ia = 0; ia < ucell.atoms[it].na; ia++) { // std::cout << "\n it = " << it << " ia = " << ia << std::endl; std::complex* sk = this->psf->get_sk(ik, it, ia, wfc_basis); //------------------------------------------------------- // Calculate G space 3D wave functions //------------------------------------------------------- - for (int iw = 0; iw < GlobalC::ucell.atoms[it].ncpp.nchi; iw++) + for (int iw = 0; iw < ucell.atoms[it].ncpp.nchi; iw++) { - if (GlobalC::ucell.atoms[it].ncpp.oc[iw] >= 0.0) + if (ucell.atoms[it].ncpp.oc[iw] >= 0.0) { - const int l = GlobalC::ucell.atoms[it].ncpp.lchi[iw]; + const int l = ucell.atoms[it].ncpp.lchi[iw]; std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, l); //----------------------------------------------------- // the factor i^l MUST BE PRESENT in order to produce @@ -324,16 +326,16 @@ void WF_atomic::atomic_wfc(const int& ik, iw, table_dimension, dq, - gk[ig].norm() * GlobalC::ucell.tpiba); + gk[ig].norm() * ucell.tpiba); } if (PARAM.inp.nspin == 4) { - if (GlobalC::ucell.atoms[it].ncpp.has_so) + if (ucell.atoms[it].ncpp.has_so) { Soc soc; soc.rot_ylm(l + 1); - const double j = GlobalC::ucell.atoms[it].ncpp.jchi[iw]; + const double j = ucell.atoms[it].ncpp.jchi[iw]; if (!(PARAM.globalv.domag || PARAM.globalv.domag_z)) { // atomic_wfc_so double fact[2]; @@ -402,10 +404,10 @@ void WF_atomic::atomic_wfc(const int& ik, } else { - for (int ib = 0; ib < GlobalC::ucell.atoms[it].ncpp.nchi; ib++) + for (int ib = 0; ib < ucell.atoms[it].ncpp.nchi; ib++) { - if ((GlobalC::ucell.atoms[it].ncpp.lchi[ib] == l) - && (fabs(GlobalC::ucell.atoms[it].ncpp.jchi[ib] - l + 0.5) < 1e-4)) + if ((ucell.atoms[it].ncpp.lchi[ib] == l) + && (fabs(ucell.atoms[it].ncpp.jchi[ib] - l + 0.5) < 1e-4)) { nc = ib; break; @@ -420,22 +422,22 @@ void WF_atomic::atomic_wfc(const int& ik, nc, table_dimension, dq, - gk[ig].norm() * GlobalC::ucell.tpiba); + gk[ig].norm() * ucell.tpiba); chiaux[ig] += flq[ig] * (l + 1.0); chiaux[ig] *= 1 / (2.0 * l + 1.0); } } // and construct the starting wavefunctions as in the noncollinear case. - // alpha = GlobalC::ucell.magnet.angle1_[it]; - // gamma = -1 * GlobalC::ucell.magnet.angle2_[it] + 0.5 * ModuleBase::PI; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamma = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; + // alpha = ucell.magnet.angle1_[it]; + // gamma = -1 * ucell.magnet.angle2_[it] + 0.5 * ModuleBase::PI; + alpha = ucell.atoms[it].angle1[ia]; + gamma = -1 * ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; for (int m = 0; m < 2 * l + 1; m++) { const int lm = l * l + m; - if (index + 2 * l + 1 > GlobalC::ucell.natomwfc) + if (index + 2 * l + 1 > ucell.natomwfc) { ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "error: too many wfcs"); } @@ -472,14 +474,14 @@ void WF_atomic::atomic_wfc(const int& ik, { // atomic_wfc_nc double alpha, gamman; std::complex fup, fdown; - // alpha = GlobalC::ucell.magnet.angle1_[it]; - // gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamman = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; + // alpha = ucell.magnet.angle1_[it]; + // gamman = -ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; + alpha = ucell.atoms[it].angle1[ia]; + gamman = -1 * ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; for (int m = 0; m < 2 * l + 1; m++) { const int lm = l * l + m; - if (index + 2 * l + 1 > GlobalC::ucell.natomwfc) + if (index + 2 * l + 1 > ucell.natomwfc) { ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "error: too many wfcs"); } @@ -545,9 +547,9 @@ void WF_atomic::atomic_wfc(const int& ik, ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wf_index", index); } - if (index != GlobalC::ucell.natomwfc) + if (index != ucell.natomwfc) { - ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "index != GlobalC::ucell.natomwfc"); + ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "index != ucell.natomwfc"); } delete[] flq; delete[] gk; diff --git a/source/module_psi/wf_atomic.h b/source/module_psi/wf_atomic.h index 2a8fb17af6..24f86cb797 100644 --- a/source/module_psi/wf_atomic.h +++ b/source/module_psi/wf_atomic.h @@ -35,16 +35,18 @@ class WF_atomic * @param sf_in [out] the structure factor * @param tab_at [out] atomic table */ - void init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at); + void init_at_1(const UnitCell& ucell, + Structure_Factor* sf_in, + ModuleBase::realArray* tab_at); - void print_PAOs() const; + void print_PAOs(const UnitCell& ucell) const; public: // template change to public, will be refactor later. added by zhengdy 20230302 int* irindex = nullptr; - void atomic_wfc(const int& ik, + void atomic_wfc(const UnitCell& ucell, + const int& ik, const int& np, - const int& lmax_wfc, const int& lmaxkb, const ModulePW::PW_Basis_K* wfc_basis, ModuleBase::ComplexMatrix& wfcatom,