From c61110069fe440e82e1bb6e487c823a1165f16ab Mon Sep 17 00:00:00 2001 From: liiutao <74701833+A-006@users.noreply.github.com> Date: Wed, 18 Dec 2024 19:29:29 +0800 Subject: [PATCH] Refactor:Remove GlobalC::ucell (#5737) * change ucell in force_stress.cpp * change ucell in force.cpp * change ucell in stress_func_cc.cpp * change ucell in rdmft_pot.cpp * change ucell in state_rmdft.cpp * change dftu_tools of ucell * remove GlobalC::ucell --- source/driver.cpp | 5 +- source/driver.h | 2 +- source/driver_run.cpp | 12 ++- source/module_esolver/esolver_ks_lcao.cpp | 4 +- .../hamilt_lcaodft/FORCE_STRESS.cpp | 12 +-- .../hamilt_lcaodft/FORCE_STRESS.h | 8 +- .../module_hamilt_lcao/module_dftu/dftu.cpp | 1 + source/module_hamilt_lcao/module_dftu/dftu.h | 1 + .../module_dftu/dftu_tools.cpp | 84 ++++++++++--------- .../hamilt_ofdft/of_stress_pw.cpp | 2 +- .../module_hamilt_pw/hamilt_pwdft/forces.cpp | 6 +- source/module_hamilt_pw/hamilt_pwdft/forces.h | 2 +- source/module_hamilt_pw/hamilt_pwdft/global.h | 1 - .../hamilt_pwdft/hamilt_pw.cpp | 2 +- .../hamilt_pwdft/stress_func.h | 3 + .../hamilt_pwdft/stress_func_cc.cpp | 47 ++++++----- .../hamilt_pwdft/stress_pw.cpp | 2 +- .../hamilt_stodft/sto_stress_pw.cpp | 4 +- .../hamilt_stodft/sto_stress_pw.h | 2 +- source/module_io/input_conv.cpp | 5 -- source/module_rdmft/rdmft.h | 4 +- source/module_rdmft/rdmft_pot.cpp | 18 ++-- source/module_rdmft/update_state_rdmft.cpp | 12 +-- 23 files changed, 129 insertions(+), 110 deletions(-) diff --git a/source/driver.cpp b/source/driver.cpp index 67476ddd5c..250ac12707 100644 --- a/source/driver.cpp +++ b/source/driver.cpp @@ -43,9 +43,6 @@ void Driver::init() // (4) close all of the running logs ModuleBase::Global_File::close_all_log(GlobalV::MY_RANK, PARAM.inp.out_alllog,PARAM.inp.calculation); - // (5) output the json file - // Json::create_Json(&GlobalC::ucell.symm,GlobalC::ucell.atoms,&INPUT); - Json::create_Json(&GlobalC::ucell, PARAM); } void Driver::print_start_info() @@ -180,7 +177,7 @@ void Driver::atomic_world() //-------------------------------------------------- // where the actual stuff is done - this->driver_run(GlobalC::ucell); + this->driver_run(); ModuleBase::timer::finish(GlobalV::ofs_running); ModuleBase::Memory::print_all(GlobalV::ofs_running); diff --git a/source/driver.h b/source/driver.h index fdcb83fade..24669e6e97 100644 --- a/source/driver.h +++ b/source/driver.h @@ -36,7 +36,7 @@ class Driver void atomic_world(); // the actual calculations - void driver_run(UnitCell& ucell); + void driver_run(); }; #endif diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 90d2c505f5..8084106fe6 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -24,13 +24,12 @@ * the configuration-changing subroutine takes force and stress and updates the * configuration */ -void Driver::driver_run(UnitCell& ucell) +void Driver::driver_run() { ModuleBase::TITLE("Driver", "driver_line"); ModuleBase::timer::tick("Driver", "driver_line"); //! 1: setup cell and atom information - // this warning should not be here, mohan 2024-05-22 #ifndef __LCAO if (PARAM.inp.basis_type == "lcao_in_pw" || PARAM.inp.basis_type == "lcao") { @@ -40,6 +39,13 @@ void Driver::driver_run(UnitCell& ucell) #endif // the life of ucell should begin here, mohan 2024-05-12 + UnitCell ucell; + ucell.setup(PARAM.inp.latname, + PARAM.inp.ntype, + PARAM.inp.lmaxmax, + PARAM.inp.init_vel, + PARAM.inp.fixed_axes); + ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); @@ -86,6 +92,8 @@ void Driver::driver_run(UnitCell& ucell) ModuleESolver::clean_esolver(p_esolver); + //! 6: output the json file + Json::create_Json(&ucell, PARAM); ModuleBase::timer::tick("Driver", "driver_line"); return; } diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 712617948a..21e2f3356e 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -302,11 +302,11 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for Force_Stress_LCAO fsl(this->RA, ucell.nat); - fsl.getForceStress(PARAM.inp.cal_force, + fsl.getForceStress(ucell, + PARAM.inp.cal_force, PARAM.inp.cal_stress, PARAM.inp.test_force, PARAM.inp.test_stress, - ucell, this->gd, this->pv, this->pelec, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 9620a5f33c..5afa1aeec7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -31,11 +31,11 @@ Force_Stress_LCAO::~Force_Stress_LCAO() { } template -void Force_Stress_LCAO::getForceStress(const bool isforce, +void Force_Stress_LCAO::getForceStress(UnitCell& ucell, + const bool isforce, const bool isstress, const bool istestf, const bool istests, - const UnitCell& ucell, const Grid_Driver& gd, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, @@ -830,7 +830,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // local pseudopotential, ewald, core correction, scc terms in force template -void Force_Stress_LCAO::calForcePwPart(const UnitCell& ucell, +void Force_Stress_LCAO::calForcePwPart(UnitCell& ucell, ModuleBase::matrix& fvl_dvl, ModuleBase::matrix& fewalds, ModuleBase::matrix& fcc, @@ -857,7 +857,7 @@ void Force_Stress_LCAO::calForcePwPart(const UnitCell& ucell, //-------------------------------------------------------- // force due to core correlation. //-------------------------------------------------------- - f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, GlobalC::ucell); + f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, ucell); //-------------------------------------------------------- // force due to self-consistent charge. //-------------------------------------------------------- @@ -975,7 +975,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn // vlocal, hartree, ewald, core correction, exchange-correlation terms in stress template -void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, +void Force_Stress_LCAO::calStressPwPart(UnitCell& ucell, ModuleBase::matrix& sigmadvl, ModuleBase::matrix& sigmahar, ModuleBase::matrix& sigmaewa, @@ -1007,7 +1007,7 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, //-------------------------------------------------------- // stress due to core correlation. //-------------------------------------------------------- - sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, nlpp.numeric, chr); + sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, nlpp.numeric, chr); //-------------------------------------------------------- // stress due to self-consistent charge. diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index c995d57b09..ecc58473a0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -29,11 +29,11 @@ class Force_Stress_LCAO Force_Stress_LCAO(Record_adj& ra, const int nat_in); ~Force_Stress_LCAO(); - void getForceStress(const bool isforce, + void getForceStress(UnitCell& ucell, + const bool isforce, const bool isstress, const bool istestf, const bool istests, - const UnitCell& ucell, const Grid_Driver& gd, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, @@ -66,7 +66,7 @@ class Force_Stress_LCAO ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm); - void calForcePwPart(const UnitCell& ucell, + void calForcePwPart(UnitCell& ucell, ModuleBase::matrix& fvl_dvl, ModuleBase::matrix& fewalds, ModuleBase::matrix& fcc, @@ -105,7 +105,7 @@ class Force_Stress_LCAO const Parallel_Orbitals& pv, const K_Vectors& kv); - void calStressPwPart(const UnitCell& ucell, + void calStressPwPart(UnitCell& ucell, ModuleBase::matrix& sigmadvl, ModuleBase::matrix& sigmahar, ModuleBase::matrix& sigmaewa, diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp index 2dd705a03c..751a9f05fd 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu.cpp @@ -59,6 +59,7 @@ void DFTU::init(UnitCell& cell, // unitcell class { orb_cutoff_ = orb->cutoffs(); } + ucell = &cell; #endif // needs reconstructions in future diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index 03a97a4905..a885340cf1 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -300,6 +300,7 @@ class DFTU void set_dmr(const elecstate::DensityMatrix, double>* dm_in_dftu_cd); private: + const UnitCell* ucell = nullptr; const elecstate::DensityMatrix* dm_in_dftu_d = nullptr; const elecstate::DensityMatrix, double>* dm_in_dftu_cd = nullptr; #endif diff --git a/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp b/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp index 363c84da89..e7d78704d7 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp @@ -12,34 +12,38 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_complex"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < this->ucell->ntype; ++it) { - if (PARAM.inp.orbital_corr[it] == -1) { + if (PARAM.inp.orbital_corr[it] == -1) + { continue; -} - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + } + for (int ia = 0; ia < this->ucell->atoms[it].na; ia++) { - const int iat = GlobalC::ucell.itia2iat(it, ia); - for (int L = 0; L <= GlobalC::ucell.atoms[it].nwl; L++) + const int iat = this->ucell->itia2iat(it, ia); + for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != PARAM.inp.orbital_corr[it]) { + if (L != PARAM.inp.orbital_corr[it]) + { continue; -} + } - for (int n = 0; n < GlobalC::ucell.atoms[it].l_nchi[L]; n++) + for (int n = 0; n < this->ucell->atoms[it].l_nchi[L]; n++) { - if (n != 0) { + if (n != 0) + { continue; -} + } for (int m1 = 0; m1 < 2 * L + 1; m1++) { 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) { + if (mu < 0) + { continue; -} + } for (int m2 = 0; m2 < 2 * L + 1; m2++) { @@ -47,13 +51,12 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); - if (nu < 0) { + if (nu < 0) + { continue; -} - + } int m1_all = m1 + (2 * L + 1) * ipol1; int m2_all = m2 + (2 * L + 1) * ipol2; - double val = get_onebody_eff_pot(it, iat, L, n, spin, m1_all, m2_all, newlocale); VU[nu * this->paraV->nrow + mu] = std::complex(val, 0.0); } // ipol2 @@ -73,43 +76,47 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_real"); ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc); - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < this->ucell->ntype; ++it) { - if (PARAM.inp.orbital_corr[it] == -1) { + if (PARAM.inp.orbital_corr[it] == -1) + { continue; -} - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + } + for (int ia = 0; ia < this->ucell->atoms[it].na; ia++) { - const int iat = GlobalC::ucell.itia2iat(it, ia); - for (int L = 0; L <= GlobalC::ucell.atoms[it].nwl; L++) + const int iat = this->ucell->itia2iat(it, ia); + for (int L = 0; L <= this->ucell->atoms[it].nwl; L++) { - if (L != PARAM.inp.orbital_corr[it]) { + if (L != PARAM.inp.orbital_corr[it]) + { continue; -} + } - for (int n = 0; n < GlobalC::ucell.atoms[it].l_nchi[L]; n++) + for (int n = 0; n < this->ucell->atoms[it].l_nchi[L]; n++) { - if (n != 0) { + if (n != 0) + { continue; -} - + } for (int m1 = 0; m1 < 2 * L + 1; m1++) { 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) { + if (mu < 0) + { continue; -} + } for (int m2 = 0; m2 < 2 * L + 1; m2++) { for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) { const int nu = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); - if (nu < 0) { + if (nu < 0) + { continue; -} + } int m1_all = m1 + (2 * L + 1) * ipol1; int m2_all = m2 + (2 * L + 1) * ipol2; @@ -157,12 +164,13 @@ double DFTU::get_onebody_eff_pot(const int T, { if (Yukawa) { - if (m0 == m1) { + if (m0 == m1) + { VU = (this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N]) * (0.5 - this->locale[iat][L][N][spin](m0, m1)); } else { VU = -(this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N]) * this->locale[iat][L][N][spin](m0, m1); -} + } } else { @@ -170,7 +178,7 @@ double DFTU::get_onebody_eff_pot(const int T, VU = (this->U[T]) * (0.5 - this->locale[iat][L][N][spin](m0, m1)); } else { VU = -(this->U[T]) * this->locale[iat][L][N][spin](m0, m1); -} + } } } else @@ -183,7 +191,7 @@ double DFTU::get_onebody_eff_pot(const int T, } else { VU = -(this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N]) * this->locale_save[iat][L][N][spin](m0, m1); -} + } } else { @@ -191,7 +199,7 @@ double DFTU::get_onebody_eff_pot(const int T, VU = (this->U[T]) * (0.5 - this->locale_save[iat][L][N][spin](m0, m1)); } else { VU = -(this->U[T]) * this->locale_save[iat][L][N][spin](m0, m1); -} + } } } diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp index d0f791d34c..c3964be8b3 100644 --- a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp @@ -79,7 +79,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, stress_loc(ucell,sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge); // nlcc - stress_cc(sigmaxcc, this->rhopw, p_sf, true, locpp.numeric, pelec->charge); + stress_cc(sigmaxcc, this->rhopw, ucell, p_sf, true, locpp.numeric, pelec->charge); // vdw term stress_vdw(sigmavdw, ucell); diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index 372836eaa9..df29f88989 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -24,7 +24,7 @@ #endif template -void Forces::cal_force(const UnitCell& ucell, +void Forces::cal_force(UnitCell& ucell, ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, @@ -161,7 +161,7 @@ void Forces::cal_force(const UnitCell& ucell, // DFT+U and DeltaSpin if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) { - this->cal_force_onsite(forceonsite, wg, wfc_basis, GlobalC::ucell, psi_in); + this->cal_force_onsite(forceonsite, wg, wfc_basis, ucell, psi_in); } } @@ -169,7 +169,7 @@ void Forces::cal_force(const UnitCell& ucell, // not relevant for PAW if (!PARAM.inp.use_paw) { - Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, GlobalC::ucell); + Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell); } else { diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h index a01aaa7f02..a300f1fba5 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.h +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.h @@ -33,7 +33,7 @@ class Forces Forces(const int nat_in) : nat(nat_in){}; ~Forces(){}; - void cal_force(const UnitCell& ucell, + void cal_force(UnitCell& ucell, ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, diff --git a/source/module_hamilt_pw/hamilt_pwdft/global.h b/source/module_hamilt_pw/hamilt_pwdft/global.h index bed170b32b..212e1e51c4 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/global.h +++ b/source/module_hamilt_pw/hamilt_pwdft/global.h @@ -265,7 +265,6 @@ namespace GlobalC #include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h" namespace GlobalC { -extern UnitCell ucell; extern Parallel_Grid Pgrid; extern Parallel_Kpoints Pkpoints; extern Restart restart; // Peize Lin add 2020.04.04 diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp index 7fe256b23d..38ccd9632c 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp @@ -118,7 +118,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, if(PARAM.inp.sc_mag_switch || PARAM.inp.dft_plus_u) { Operator* onsite_proj - = new OnsiteProj>(isk, &GlobalC::ucell, PARAM.inp.sc_mag_switch, (PARAM.inp.dft_plus_u>0)); + = new OnsiteProj>(isk, ucell, PARAM.inp.sc_mag_switch, (PARAM.inp.dft_plus_u>0)); this->ops->add(onsite_proj); } return; diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h index a81dbc9d93..878206ad38 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h @@ -114,12 +114,15 @@ class Stress_Func // 5) the stress from the non-linear core correction (if any) void stress_cc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, + UnitCell& ucell, const Structure_Factor* p_sf, const bool is_pw, const bool *numeric, const Charge* const chr); // nonlinear core correction stress in PW or LCAO basis void deriv_drhoc(const bool& numeric, + const double& omega, + const double& tpiba2, const int mesh, const FPTYPE* r, const FPTYPE* rab, diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index 63f892c0c2..ab8d9b3fa1 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -15,6 +15,7 @@ template void Stress_Func::stress_cc(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, + UnitCell& ucell, const Structure_Factor* p_sf, const bool is_pw, const bool *numeric, @@ -34,9 +35,9 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, FPTYPE* rhocg; int judge=0; - for(int nt=0;nt::stress_cc(ModuleBase::matrix& sigma, { #ifdef USE_LIBXC const auto etxc_vtxc_v - = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, GlobalC::ucell.omega, GlobalC::ucell.tpiba, chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell.omega, ucell.tpiba, chr); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); @@ -65,8 +66,8 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, } else { - elecstate::cal_ux(GlobalC::ucell); - const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &GlobalC::ucell); + elecstate::cal_ux(ucell); + const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell); // etxc = std::get<0>(etxc_vtxc_v); // may delete? // vtxc = std::get<1>(etxc_vtxc_v); // may delete? vxc = std::get<2>(etxc_vtxc_v); @@ -103,17 +104,19 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, rhocg= new FPTYPE [rho_basis->ngg]; sigmadiag=0.0; - for(int nt=0;ntderiv_drhoc( numeric, - GlobalC::ucell.atoms[nt].ncpp.msh, - GlobalC::ucell.atoms[nt].ncpp.r.data(), - GlobalC::ucell.atoms[nt].ncpp.rab.data(), - GlobalC::ucell.atoms[nt].ncpp.rho_atc.data(), + ucell.omega, + ucell.tpiba2, + ucell.atoms[nt].ncpp.msh, + ucell.atoms[nt].ncpp.r.data(), + ucell.atoms[nt].ncpp.rab.data(), + ucell.atoms[nt].ncpp.rho_atc.data(), rhocg, rho_basis, 1); @@ -135,10 +138,12 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, } this->deriv_drhoc ( numeric, - GlobalC::ucell.atoms[nt].ncpp.msh, - GlobalC::ucell.atoms[nt].ncpp.r.data(), - GlobalC::ucell.atoms[nt].ncpp.rab.data(), - GlobalC::ucell.atoms[nt].ncpp.rho_atc.data(), + ucell.omega, + ucell.tpiba2, + ucell.atoms[nt].ncpp.msh, + ucell.atoms[nt].ncpp.r.data(), + ucell.atoms[nt].ncpp.rab.data(), + ucell.atoms[nt].ncpp.rho_atc.data(), rhocg, rho_basis, 0); @@ -162,7 +167,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, { const std::complex t = conj(psic[ig]) * p_sf->strucFac(nt, ig) * rhocg[rho_basis->ig2igg[ig]] - * GlobalC::ucell.tpiba * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / norm_g * fact; + * ucell.tpiba * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / norm_g * fact; // sigmacc [l][ m] += t.real(); local_sigma(l,m) += t.real(); }//end m @@ -209,6 +214,8 @@ template void Stress_Func::deriv_drhoc ( const bool &numeric, + const double& omega, + const double& tpiba2, const int mesh, const FPTYPE *r, const FPTYPE *rab, @@ -254,7 +261,7 @@ void Stress_Func::deriv_drhoc aux [ir] = r [ir] * r [ir] * rhoc [ir]; } ModuleBase::Integral::Simpson_Integral(mesh, aux.data(), rab, rhocg1); - drhocg [0] = ModuleBase::FOUR_PI * rhocg1 / GlobalC::ucell.omega; + drhocg [0] = ModuleBase::FOUR_PI * rhocg1 / omega; igl0 = 1; } else @@ -273,7 +280,7 @@ void Stress_Func::deriv_drhoc #endif for(int igl = igl0;igl< rho_basis->ngg;igl++) { - gx_arr[igl] = sqrt(rho_basis->gg_uniq[igl] * GlobalC::ucell.tpiba2); + gx_arr[igl] = sqrt(rho_basis->gg_uniq[igl] * tpiba2); } double *r_d = nullptr; @@ -298,12 +305,12 @@ void Stress_Func::deriv_drhoc if(this->device == base_device::GpuDevice) { hamilt::cal_stress_drhoc_aux_op()( - r_d,rhoc_d,gx_arr_d+igl0,rab_d,drhocg_d+igl0,mesh,igl0,rho_basis->ngg-igl0,GlobalC::ucell.omega,type); + r_d,rhoc_d,gx_arr_d+igl0,rab_d,drhocg_d+igl0,mesh,igl0,rho_basis->ngg-igl0,omega,type); syncmem_var_d2h_op()(this->cpu_ctx, this->ctx, drhocg+igl0, drhocg_d+igl0, rho_basis->ngg-igl0); } else { hamilt::cal_stress_drhoc_aux_op()( - r,rhoc,gx_arr.data()+igl0,rab,drhocg+igl0,mesh,igl0,rho_basis->ngg-igl0,GlobalC::ucell.omega,type); + r,rhoc,gx_arr.data()+igl0,rab,drhocg+igl0,mesh,igl0,rho_basis->ngg-igl0,omega,type); } delmem_var_op()(this->ctx, r_d); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp index e9cc5ded2b..33bf039a62 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp @@ -97,7 +97,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->stress_loc(ucell,sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge); // nlcc - this->stress_cc(sigmaxcc, rho_basis, p_sf, 1, nlpp.numeric, pelec->charge); + this->stress_cc(sigmaxcc, rho_basis, ucell, p_sf, 1, nlpp.numeric, pelec->charge); // nonlocal this->stress_nl(sigmanl, this->pelec->wg, this->pelec->ekb, p_sf, p_kv, p_symm, wfc_basis, d_psi_in, nlpp, ucell); 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 f44172c52a..9115b5a053 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp @@ -20,7 +20,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, const Stochastic_WF, Device>& stowf, const Charge* const chr, pseudopot_cell_vnl* nlpp, - const UnitCell& ucell_in) + UnitCell& ucell_in) { ModuleBase::TITLE("Sto_Stress_PW", "cal_stress"); ModuleBase::timer::tick("Sto_Stress_PW", "cal_stress"); @@ -55,7 +55,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->stress_loc(ucell_in,sigmaloc, rho_basis, nlpp->vloc, p_sf, true, chr); // nlcc - this->stress_cc(sigmaxcc, rho_basis, p_sf, true, nlpp->numeric, chr); + this->stress_cc(sigmaxcc, rho_basis, ucell_in, p_sf, true, nlpp->numeric, chr); // nonlocal this->sto_stress_nl(sigmanl, wg, p_sf, p_symm, p_kv, wfc_basis, *nlpp, ucell_in, psi_in, stowf); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h index 7b25166dac..4e67d74110 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h +++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h @@ -26,7 +26,7 @@ class Sto_Stress_PW : public Stress_Func const Stochastic_WF, Device>& stowf, const Charge* const chr, pseudopot_cell_vnl* nlpp_in, - const UnitCell& ucell_in); + UnitCell& ucell_in); private: void sto_stress_kin(ModuleBase::matrix& sigma, diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 7ce7f0d764..65ff27dbc1 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -170,11 +170,6 @@ void Input_Conv::Convert() ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "pseudo_dir", PARAM.inp.pseudo_dir); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "orbital_dir", PARAM.inp.orbital_dir); // GlobalV::global_pseudo_type = PARAM.inp.pseudo_type; - GlobalC::ucell.setup(PARAM.inp.latname, - PARAM.inp.ntype, - PARAM.inp.lmaxmax, - PARAM.inp.init_vel, - PARAM.inp.fixed_axes); if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") { diff --git a/source/module_rdmft/rdmft.h b/source/module_rdmft/rdmft.h index 423fc7b94a..bc02480637 100644 --- a/source/module_rdmft/rdmft.h +++ b/source/module_rdmft/rdmft.h @@ -99,7 +99,7 @@ class RDMFT //! update in elec-step // Or we can use rdmft_solver.wfc/occ_number directly when optimizing, so that the update_elec() function does not require parameters. - void update_elec(const UnitCell& ucell, const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in = nullptr); + void update_elec(UnitCell& ucell, const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in = nullptr); //! obtain the gradient of total energy with respect to occupation number and wfc double cal_E_grad_wfc_occ_num(); @@ -131,7 +131,7 @@ class RDMFT //! get the total Hamilton in k-space void cal_Hk_Hpsi(); - void update_charge(); + void update_charge(UnitCell& ucell); private: diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index 06d9c5a3ff..e576c17942 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -54,7 +54,7 @@ void RDMFT::cal_V_TV() V_ekinetic_potential = new hamilt::EkineticNew>(hsk_TV, kv->kvec_d, HR_TV, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, two_center_bundle->kinetic_orb.get()); @@ -62,7 +62,7 @@ void RDMFT::cal_V_TV() V_nonlocal = new hamilt::NonlocalNew>(hsk_TV, kv->kvec_d, HR_TV, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, two_center_bundle->overlap_orb_beta.get()); @@ -74,7 +74,7 @@ void RDMFT::cal_V_TV() kv->kvec_d, this->pelec->pot, HR_TV, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, @@ -91,7 +91,7 @@ void RDMFT::cal_V_TV() kv->kvec_d, this->pelec->pot, HR_TV, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, @@ -122,7 +122,7 @@ void RDMFT::cal_V_hartree() kv->kvec_d, this->pelec->pot, HR_hartree, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, @@ -140,7 +140,7 @@ void RDMFT::cal_V_hartree() kv->kvec_d, this->pelec->pot, HR_hartree, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, @@ -170,7 +170,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) // elecstate::DensityMatrix DM_test(ParaV, nspin, kv->kvec_d, nk_total); // elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test); - // DM_test.init_DMR(this->gd, &GlobalC::ucell); + // DM_test.init_DMR(this->gd, this->ucell); // DM_test.cal_DMR(); // // compare DM_XC and DM get in update_charge(or ABACUS) @@ -202,7 +202,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) kv->kvec_d, this->pelec->pot, HR_dft_XC, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, @@ -222,7 +222,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) kv->kvec_d, this->pelec->pot, HR_dft_XC, - &GlobalC::ucell, + this->ucell, orb->cutoffs(), this->gd, nspin, diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index db5f47b6da..abe56d71c3 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -47,7 +47,7 @@ void RDMFT::update_ion(UnitCell& ucell_in, template -void RDMFT::update_elec(const UnitCell& ucell, +void RDMFT::update_elec(UnitCell& ucell, const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in) { @@ -71,7 +71,7 @@ void RDMFT::update_elec(const UnitCell& ucell, } // update charge - this->update_charge(); + this->update_charge(ucell); // "default" = "pbe" // if( !only_exx_type || this->cal_E_type != 1 ) @@ -91,14 +91,14 @@ void RDMFT::update_elec(const UnitCell& ucell, // this code is copying from function ElecStateLCAO::psiToRho(), in elecstate_lcao.cpp template -void RDMFT::update_charge() +void RDMFT::update_charge(UnitCell& ucell) { if( PARAM.inp.gamma_only ) { // calculate DMK and DMR elecstate::DensityMatrix DM_gamma_only(ParaV, nspin); elecstate::cal_dm_psi(ParaV, wg, wfc, DM_gamma_only); - DM_gamma_only.init_DMR(this->gd, &GlobalC::ucell); + DM_gamma_only.init_DMR(this->gd, &ucell); DM_gamma_only.cal_DMR(); for (int is = 0; is < nspin; is++) @@ -128,7 +128,7 @@ void RDMFT::update_charge() // calculate DMK and DMR elecstate::DensityMatrix DM(ParaV, nspin, kv->kvec_d, nk_total); elecstate::cal_dm_psi(ParaV, wg, wfc, DM); - DM.init_DMR(this->gd, &GlobalC::ucell); + DM.init_DMR(this->gd, &ucell); DM.cal_DMR(); for (int is = 0; is < nspin; is++) @@ -160,7 +160,7 @@ void RDMFT::update_charge() Symmetry_rho srho; for (int is = 0; is < nspin; is++) { - srho.begin(is, *(this->charge), rho_basis, GlobalC::ucell.symm); + srho.begin(is, *(this->charge), rho_basis, ucell.symm); } }