diff --git a/source/module_cell/setup_nonlocal.cpp b/source/module_cell/setup_nonlocal.cpp index 9fb7280c1f..f62982db49 100644 --- a/source/module_cell/setup_nonlocal.cpp +++ b/source/module_cell/setup_nonlocal.cpp @@ -34,7 +34,7 @@ void InfoNonlocal::Set_NonLocal(const int& it, ModuleBase::TITLE("InfoNonlocal", "Set_NonLocal"); // set a pointer - // Atom* atom = &GlobalC::ucell.atoms[it]; + // Atom* atom = &ucell.atoms[it]; // get the number of non-local projectors n_projectors = atom->ncpp.nbeta; diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index c84ce504d2..64b0db0ae9 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -207,6 +207,7 @@ void ElecState::calEBand() void ElecState::init_scf(const int istep, + const UnitCell& ucell, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, ModuleSymmetry::Symmetry& symm, @@ -215,7 +216,7 @@ void ElecState::init_scf(const int istep, //! core correction potential. if (!PARAM.inp.use_paw) { - this->charge->set_rho_core(strucfac, numeric); + this->charge->set_rho_core(ucell,strucfac, numeric); } else { @@ -226,7 +227,7 @@ void ElecState::init_scf(const int istep, // choose charge density from ionic step 0. if (istep == 0) { - this->charge->init_rho(this->eferm, strucfac, symm, (const void*)this->klist, wfcpw); + this->charge->init_rho(ucell,this->eferm, strucfac, symm, (const void*)this->klist, wfcpw); this->charge->check_rho(); // check the rho } diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index a90555a249..7bb874af60 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -109,6 +109,7 @@ class ElecState * @param wfcpw PW basis for wave function if needed */ void init_scf(const int istep, + const UnitCell& ucell, const ModuleBase::ComplexMatrix& strucfac, const bool* numeric, ModuleSymmetry::Symmetry& symm, @@ -124,9 +125,9 @@ class ElecState public: // something aboud energies. See elecstate_energy.cpp void cal_bandgap(); - void cal_bandgap_updw(); + void cal_bandgap_updw(const UnitCell& ucell); - double cal_delta_eband() const; + double cal_delta_eband(const UnitCell& ucell) const; double cal_delta_escf() const; ModuleBase::matrix vnew; @@ -173,7 +174,8 @@ class ElecState ModuleBase::matrix wg; ///< occupation weight for each k-point and band public: // print something. See elecstate_print.cpp - void print_etot(const bool converged, + void print_etot(const UnitCell& ucell, + const bool converged, const int& iter, const double& scf_thr, const double& scf_thr_kin, diff --git a/source/module_elecstate/elecstate_energy.cpp b/source/module_elecstate/elecstate_energy.cpp index 86a02d7364..32f8bc9b2d 100644 --- a/source/module_elecstate/elecstate_energy.cpp +++ b/source/module_elecstate/elecstate_energy.cpp @@ -43,7 +43,7 @@ void ElecState::cal_bandgap() /// @brief calculate spin up & down band gap /// @todo add isk[ik] so as to discriminate different spins -void ElecState::cal_bandgap_updw() +void ElecState::cal_bandgap_updw(const UnitCell& ucell) { if (this->ekb.nr == 0 || this->ekb.nc == 0) { // which means no homo and no lumo @@ -90,7 +90,7 @@ void ElecState::cal_bandgap_updw() } /// @brief calculate deband -double ElecState::cal_delta_eband() const +double ElecState::cal_delta_eband(const UnitCell& ucell) const { // out potentials from potential mixing // total energy and band energy corrections @@ -109,7 +109,7 @@ double ElecState::cal_delta_eband() const { ModuleBase::matrix v_xc; const std::tuple etxc_vtxc_v - = XC_Functional::v_xc(this->charge->nrxx, this->charge, &GlobalC::ucell); + = XC_Functional::v_xc(this->charge->nrxx, this->charge, &ucell); v_xc = std::get<2>(etxc_vtxc_v); for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) diff --git a/source/module_elecstate/elecstate_getters.cpp b/source/module_elecstate/elecstate_getters.cpp index 30bf81f94c..d39cf7ec42 100644 --- a/source/module_elecstate/elecstate_getters.cpp +++ b/source/module_elecstate/elecstate_getters.cpp @@ -8,54 +8,13 @@ namespace elecstate { -double get_ucell_omega() -{ - return GlobalC::ucell.omega; -} - -double get_ucell_tpiba() -{ - return GlobalC::ucell.tpiba; -} int get_xc_func_type() { return XC_Functional::get_func_type(); } -std::string get_input_vdw_method() -{ - return PARAM.inp.vdw_method; -} - -double get_ucell_tot_magnetization() -{ - return GlobalC::ucell.magnet.tot_magnetization; -} -double get_ucell_abs_magnetization() -{ - return GlobalC::ucell.magnet.abs_magnetization; -} -double get_ucell_tot_magnetization_nc_x() -{ - return GlobalC::ucell.magnet.tot_magnetization_nc[0]; -} - -double get_ucell_tot_magnetization_nc_y() -{ - return GlobalC::ucell.magnet.tot_magnetization_nc[1]; -} - -double get_ucell_tot_magnetization_nc_z() -{ - return GlobalC::ucell.magnet.tot_magnetization_nc[2]; -} - -std::string get_ks_solver_type() -{ - return PARAM.inp.ks_solver; -} } // namespace elecstate diff --git a/source/module_elecstate/elecstate_getters.h b/source/module_elecstate/elecstate_getters.h index 4ee4ca57c0..6c412b4e00 100644 --- a/source/module_elecstate/elecstate_getters.h +++ b/source/module_elecstate/elecstate_getters.h @@ -7,26 +7,8 @@ namespace elecstate { -/// @brief get the value of GlobalC::ucell.omega -double get_ucell_omega(); -/// @brief get the value of GlobalC::ucell.tpiba -double get_ucell_tpiba(); /// @brief get the value of XC_Functional::func_type int get_xc_func_type(); -/// @brief get the value of INPUT.vdw_method -std::string get_input_vdw_method(); -/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization -double get_ucell_tot_magnetization(); -/// @brief get the value of GlobalC::ucell.magnet.abs_magnetization -double get_ucell_abs_magnetization(); -/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[0] -double get_ucell_tot_magnetization_nc_x(); -/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[1] -double get_ucell_tot_magnetization_nc_y(); -/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[2] -double get_ucell_tot_magnetization_nc_z(); -/// @brief get the type of KS_SOLVER -std::string get_ks_solver_type(); } // namespace elecstate diff --git a/source/module_elecstate/elecstate_print.cpp b/source/module_elecstate/elecstate_print.cpp index 2d738345b3..3af86843ef 100644 --- a/source/module_elecstate/elecstate_print.cpp +++ b/source/module_elecstate/elecstate_print.cpp @@ -289,7 +289,8 @@ void ElecState::print_band(const int& ik, const int& printe, const int& iter) /// @param pw_diag_thr: threshold for diagonalization /// @param avg_iter: averaged diagonalization iteration of each scf iteration /// @param print: if print to screen -void ElecState::print_etot(const bool converged, +void ElecState::print_etot(const UnitCell& ucell, + const bool converged, const int& iter_in, const double& scf_thr, const double& scf_thr_kin, @@ -339,7 +340,7 @@ void ElecState::print_etot(const bool converged, energies_Ry.push_back(this->f_en.demet); titles.push_back("E_descf"); energies_Ry.push_back(this->f_en.descf); - std::string vdw_method = get_input_vdw_method(); + std::string vdw_method = PARAM.inp.vdw_method; if (vdw_method == "d2") // Peize Lin add 2014-04, update 2021-03-09 { titles.push_back("E_vdwD2"); @@ -429,13 +430,13 @@ void ElecState::print_etot(const bool converged, switch (PARAM.inp.nspin) { case 2: - mag = {get_ucell_tot_magnetization(), get_ucell_abs_magnetization()}; + mag = {ucell.magnet.tot_magnetization, ucell.magnet.abs_magnetization}; break; case 4: - mag = {get_ucell_tot_magnetization_nc_x(), - get_ucell_tot_magnetization_nc_y(), - get_ucell_tot_magnetization_nc_z(), - get_ucell_abs_magnetization()}; + mag = {ucell.magnet.tot_magnetization_nc[0], + ucell.magnet.tot_magnetization_nc[1], + ucell.magnet.tot_magnetization_nc[2], + ucell.magnet.abs_magnetization}; break; default: mag = {}; @@ -446,7 +447,7 @@ void ElecState::print_etot(const bool converged, { drho.push_back(scf_thr_kin); } - elecstate::print_scf_iterinfo(get_ks_solver_type(), + elecstate::print_scf_iterinfo(PARAM.inp.ks_solver, iter, 6, mag, diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index 696a3a75a2..05fb2b57f3 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -176,7 +176,7 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) this->basis->recip_to_real(this->ctx, &psi(ibnd,npwx), this->wfcr_another_spin, ik); - const auto w1 = static_cast(this->wg(ik, ibnd) / get_ucell_omega()); + const auto w1 = static_cast(this->wg(ik, ibnd) / ucell->omega); if (w1 != 0.0) { @@ -202,7 +202,7 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) this->basis->recip_to_real(this->ctx, &psi(ibnd,0), this->wfcr, ik); - const auto w1 = static_cast(this->wg(ik, ibnd) / get_ucell_omega()); + const auto w1 = static_cast(this->wg(ik, ibnd) / ucell->omega); if (w1 != 0.0) { @@ -222,7 +222,7 @@ void ElecStatePW::rhoBandK(const psi::Psi& psi) j, npw, this->basis->npwk_max, - static_cast(get_ucell_tpiba()), + static_cast(ucell->tpiba), this->basis->template get_gcar_data(), this->basis->template get_kvec_c_data(), &psi(ibnd, 0), diff --git a/source/module_elecstate/elecstate_pw_cal_tau.cpp b/source/module_elecstate/elecstate_pw_cal_tau.cpp index 0662bb4425..fd07f834af 100644 --- a/source/module_elecstate/elecstate_pw_cal_tau.cpp +++ b/source/module_elecstate/elecstate_pw_cal_tau.cpp @@ -26,7 +26,7 @@ void ElecStatePW::cal_tau(const psi::Psi& psi) { this->basis->recip_to_real(this->ctx, &psi(ibnd,0), this->wfcr, ik); - const auto w1 = static_cast(this->wg(ik, ibnd) / get_ucell_omega()); + const auto w1 = static_cast(this->wg(ik, ibnd) / ucell->omega); // kinetic energy density for (int j = 0; j < 3; j++) @@ -38,7 +38,7 @@ void ElecStatePW::cal_tau(const psi::Psi& psi) j, npw, this->basis->npwk_max, - static_cast(get_ucell_tpiba()), + static_cast(ucell->tpiba), this->basis->template get_gcar_data(), this->basis->template get_kvec_c_data(), &psi(ibnd, 0), diff --git a/source/module_elecstate/magnetism.cpp b/source/module_elecstate/magnetism.cpp index 89bf71eb38..7b598043f0 100644 --- a/source/module_elecstate/magnetism.cpp +++ b/source/module_elecstate/magnetism.cpp @@ -15,7 +15,11 @@ Magnetism::~Magnetism() delete[] this->start_magnetization; } -void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const double* const * rho, double* nelec_spin) +void Magnetism::compute_magnetization(const double omega, + const int& nrxx, + const int& nxyz, + const double* const * rho, + double* nelec_spin) { if (PARAM.inp.nspin==2) { @@ -32,8 +36,8 @@ void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const do Parallel_Reduce::reduce_pool(this->tot_magnetization); Parallel_Reduce::reduce_pool(this->abs_magnetization); #endif - this->tot_magnetization *= elecstate::get_ucell_omega() / nxyz; - this->abs_magnetization *= elecstate::get_ucell_omega() / nxyz; + this->tot_magnetization *= omega / nxyz; + this->abs_magnetization *= omega / nxyz; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"total magnetism (Bohr mag/cell)",this->tot_magnetization); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization); @@ -65,8 +69,8 @@ void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const do Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3); Parallel_Reduce::reduce_pool(this->abs_magnetization); #endif - for(int i=0;i<3;i++)this->tot_magnetization_nc[i] *= elecstate::get_ucell_omega() / nxyz; - this->abs_magnetization *= elecstate::get_ucell_omega() / nxyz; + for(int i=0;i<3;i++)this->tot_magnetization_nc[i] *= omega / nxyz; + this->abs_magnetization *= omega / nxyz; GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<tot_magnetization_nc[0]<<'\t'<tot_magnetization_nc[1]<<'\t'<tot_magnetization_nc[2]<<'\n'; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization); } diff --git a/source/module_elecstate/magnetism.h b/source/module_elecstate/magnetism.h index fbc750b10c..3d1c8c1859 100644 --- a/source/module_elecstate/magnetism.h +++ b/source/module_elecstate/magnetism.h @@ -21,7 +21,7 @@ class Magnetism double tot_magnetization_nc[3]; double abs_magnetization; - void compute_magnetization(const int& nrxx, const int& nxyz, const double* const * rho, double* nelec_spin = nullptr); + void compute_magnetization(const double omega, const int& nrxx, const int& nxyz, const double* const * rho, double* nelec_spin = nullptr); ModuleBase::Vector3 *m_loc_; //magnetization for each element along c-axis double *angle1_; //angle between c-axis and real spin std::vector diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index 6a2405b579..a15143ce6e 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -218,7 +218,7 @@ double Charge::sum_rho() const } // multiply the sum of charge density by a factor - sum_rho *= elecstate::get_ucell_omega() / static_cast(this->rhopw->nxyz); + sum_rho *= this->omega / static_cast(this->rhopw->nxyz); #ifdef __MPI Parallel_Reduce::reduce_pool(sum_rho); @@ -722,7 +722,7 @@ double Charge::cal_rho2ne(const double* rho_in) const #ifdef __MPI Parallel_Reduce::reduce_pool(ne); #endif - ne = ne * elecstate::get_ucell_omega() / (double)this->rhopw->nxyz; + ne = ne * this->omega / (double)this->rhopw->nxyz; return ne; } diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index bc8a1aa4ec..6b214b716c 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -70,7 +70,8 @@ class Charge * @param klist [in] k points list if needed * @param wfcpw [in] PW basis for wave function if needed */ - void init_rho(elecstate::efermi& eferm_iout, + void init_rho(const UnitCell& ucell, + elecstate::efermi& eferm_iout, const ModuleBase::ComplexMatrix& strucFac, ModuleSymmetry::Symmetry& symm, const void* klist = nullptr, @@ -84,7 +85,9 @@ class Charge const ModuleBase::ComplexMatrix& strucFac, const UnitCell& ucell) const; - void set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric); + void set_rho_core(const UnitCell& ucell, + const ModuleBase::ComplexMatrix& structure_factor, + const bool* numeric); void set_rho_core_paw(); void renormalize_rho(); @@ -97,6 +100,8 @@ class Charge void non_linear_core_correction ( const bool &numeric, + const double omega, + const double tpiba2, const int mesh, const double *r, const double *rab, @@ -144,7 +149,8 @@ class Charge void destroy(); // free arrays liuyu 2023-03-12 bool allocate_rho; - + double omega; // volume of the unit cell + double tpiba2; // 2*pi/lat0 bool allocate_rho_final_scf; // LiuXh add 20180606 #ifdef __MPI private: diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp index b38c5cdd5c..2f6953bf9c 100644 --- a/source/module_elecstate/module_charge/charge_init.cpp +++ b/source/module_elecstate/module_charge/charge_init.cpp @@ -21,7 +21,8 @@ #include "module_cell/module_paw/paw_cell.h" #endif -void Charge::init_rho(elecstate::efermi& eferm_iout, +void Charge::init_rho(const UnitCell& ucell, + elecstate::efermi& eferm_iout, const ModuleBase::ComplexMatrix& strucFac, ModuleSymmetry::Symmetry& symm, const void* klist, @@ -31,6 +32,8 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, std::cout << " START CHARGE : " << PARAM.inp.init_chg << std::endl; bool read_error = false; + this->omega = ucell.omega; + this->tpiba2 = ucell.tpiba2; if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto") { GlobalV::ofs_running << " try to read charge from file : " << std::endl; @@ -58,7 +61,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, GlobalV::ofs_running, ssc.str(), this->rho[is], - GlobalC::ucell.nat)) + ucell.nat)) { GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl; } @@ -108,7 +111,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, GlobalV::ofs_running, ssc.str(), this->kin_r[is], - GlobalC::ucell.nat)) + ucell.nat)) { GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl; } @@ -134,7 +137,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, if (PARAM.inp.init_chg == "atomic" || (PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17 { - this->atomic_rho(PARAM.inp.nspin, GlobalC::ucell.omega, rho, strucFac, GlobalC::ucell); + this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell); // liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi @@ -171,7 +174,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, GlobalV::ofs_running, ssc.str(), this->rho[is], - GlobalC::ucell.nat)) + ucell.nat)) { GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl; } @@ -199,7 +202,9 @@ void Charge::init_rho(elecstate::efermi& eferm_iout, //========================================================== // computes the core charge on the real space 3D mesh. //========================================================== -void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric) +void Charge::set_rho_core(const UnitCell& ucell, + const ModuleBase::ComplexMatrix& structure_factor, + const bool* numeric) { ModuleBase::TITLE("Charge","set_rho_core"); ModuleBase::timer::tick("Charge","set_rho_core"); @@ -216,9 +221,9 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con // int ig = 0; bool bl = false; - for (int it = 0; it *vg = new std::complex[this->rhopw->npw]; - for (int it = 0; it < GlobalC::ucell.ntype;it++) + for (int it = 0; it < ucell.ntype;it++) { - if (GlobalC::ucell.atoms[it].ncpp.nlcc) + if (ucell.atoms[it].ncpp.nlcc) { //---------------------------------------------------------- // EXPLAIN : drhoc compute the radial fourier transform for @@ -248,10 +253,12 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con //---------------------------------------------------------- this->non_linear_core_correction( numeric, - GlobalC::ucell.atoms[it].ncpp.msh, - GlobalC::ucell.atoms[it].ncpp.r.data(), - GlobalC::ucell.atoms[it].ncpp.rab.data(), - GlobalC::ucell.atoms[it].ncpp.rho_atc.data(), + ucell.omega, + ucell.tpiba2, + ucell.atoms[it].ncpp.msh, + ucell.atoms[it].ncpp.r.data(), + ucell.atoms[it].ncpp.rab.data(), + ucell.atoms[it].ncpp.rho_atc.data(), rhocg); //---------------------------------------------------------- // EXPLAIN : multiply by the structure factor and sum @@ -296,8 +303,8 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con // mohan changed 2010-2-2, make this same as in atomic_rho. // still lack something...... - rhoneg /= this->rhopw->nxyz * GlobalC::ucell.omega; - rhoima /= this->rhopw->nxyz * GlobalC::ucell.omega; + rhoneg /= this->rhopw->nxyz * ucell.omega; + rhoima /= this->rhopw->nxyz * ucell.omega; // calculate core_only exch-corr energy etxcc=E_xc[rho_core] if required // The term was present in previous versions of the code but it shouldn't @@ -322,6 +329,8 @@ void Charge::set_rho_core_paw() void Charge::non_linear_core_correction ( const bool &numeric, + const double omega, + const double tpiba2, const int mesh, const double *r, const double *rab, @@ -356,7 +365,7 @@ void Charge::non_linear_core_correction } ModuleBase::Integral::Simpson_Integral(mesh, aux, rab, rhocg1); //rhocg [1] = fpi * rhocg1 / omega; - rhocg [0] = ModuleBase::FOUR_PI * rhocg1 / GlobalC::ucell.omega;//mohan modify 2008-01-19 + rhocg [0] = ModuleBase::FOUR_PI * rhocg1 / omega;//mohan modify 2008-01-19 } igl0 = 1; } @@ -370,14 +379,14 @@ void Charge::non_linear_core_correction // G <> 0 term for (int igl = igl_beg; igl < igl_end;igl++) { - gx = sqrt(this->rhopw->gg_uniq[igl] * GlobalC::ucell.tpiba2); + gx = sqrt(this->rhopw->gg_uniq[igl] * tpiba2); ModuleBase::Sphbes::Spherical_Bessel(mesh, r, gx, 0, aux); for (int ir = 0;ir < mesh; ir++) { aux [ir] = r[ir] * r[ir] * rhoc [ir] * aux [ir]; } // enddo ModuleBase::Integral::Simpson_Integral(mesh, aux, rab, rhocg1); - rhocg [igl] = ModuleBase::FOUR_PI * rhocg1 / GlobalC::ucell.omega; + rhocg [igl] = ModuleBase::FOUR_PI * rhocg1 / omega; } // enddo delete [] aux; } diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp index 37860691b7..81960fda3c 100644 --- a/source/module_elecstate/module_charge/charge_mixing.cpp +++ b/source/module_elecstate/module_charge/charge_mixing.cpp @@ -25,7 +25,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, const double& mixing_gg0_mag_in, const double& mixing_gg0_min_in, const double& mixing_angle_in, - const bool& mixing_dmr_in) + const bool& mixing_dmr_in, + const double& omega_in, + const double& tpiba_in) { // get private mixing parameters this->mixing_mode = mixing_mode_in; @@ -38,7 +40,10 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, this->mixing_gg0_min = mixing_gg0_min_in; this->mixing_angle = mixing_angle_in; this->mixing_dmr = mixing_dmr_in; - + this->omega = omega_in; + this->tpiba = tpiba_in; + this->tpiba2 = tpiba_in * tpiba_in; + // check the paramters if (this->mixing_beta > 1.0 || this->mixing_beta < 0.0) { diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index 9e85908d12..d1e56dd02b 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -43,7 +43,9 @@ class Charge_Mixing const double& mixing_gg0_mag_in, const double& mixing_gg0_min_in, const double& mixing_angle_in, - const bool& mixing_dmr_in); + const bool& mixing_dmr_in, + const double& omega_in, + const double& tpiba_in); void close_kerker_gg0() { mixing_gg0 = 0.0; mixing_gg0_mag = 0.0; } /** @@ -75,8 +77,8 @@ class Charge_Mixing * @brief Get the drho between rho and rho_save, similar for get_dkin * */ - double get_drho(Charge* chr, const double nelec); - double get_dkin(Charge* chr, const double nelec); + double get_drho(Charge* chr, const double nelec , const double omega); + double get_dkin(Charge* chr, const double nelec , const double omega); /** * @brief reset mixing, actually we only call init_mixing() to reset mixing instead of this function @@ -129,6 +131,9 @@ class Charge_Mixing double mixing_gg0_min = 0.1; ///< minimum kerker coefficient double mixing_angle = 0.0; ///< mixing angle for nspin=4 bool mixing_dmr = false; ///< whether to mixing real space density matrix + double omega = 0.0; ///< omega for mixing + double tpiba = 0.0; ///< 2*pi/lat0 + double tpiba2 = 0.0; ///< 2*pi/lat0^2 std::vector _drho_history; ///< history of drho used to determine the oscillation, size is scf_nmax diff --git a/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp b/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp index 3564b86f98..c0999c824b 100644 --- a/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp +++ b/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp @@ -45,7 +45,7 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) amin = this->mixing_beta; } - gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); + gg0 = std::pow(fac * 0.529177 / tpiba2, 2); #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif @@ -110,7 +110,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) amin = this->mixing_beta; } - gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2); + gg0 = std::pow(fac * 0.529177 / tpiba, 2); #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif diff --git a/source/module_elecstate/module_charge/charge_mixing_residual.cpp b/source/module_elecstate/module_charge/charge_mixing_residual.cpp index d2683629d5..9597eae558 100644 --- a/source/module_elecstate/module_charge/charge_mixing_residual.cpp +++ b/source/module_elecstate/module_charge/charge_mixing_residual.cpp @@ -5,7 +5,9 @@ #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_base/parallel_reduce.h" -double Charge_Mixing::get_drho(Charge* chr, const double nelec) +double Charge_Mixing::get_drho(Charge* chr, + const double nelec, + const double omega) { ModuleBase::TITLE("Charge_Mixing", "get_drho"); ModuleBase::timer::tick("Charge_Mixing", "get_drho"); @@ -60,9 +62,9 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) Parallel_Reduce::reduce_pool(drho); #endif assert(nelec != 0); - assert(GlobalC::ucell.omega > 0); + assert(omega > 0); assert(this->rhopw->nxyz > 0); - drho *= GlobalC::ucell.omega / static_cast(this->rhopw->nxyz); + drho *= omega / static_cast(this->rhopw->nxyz); drho /= nelec; } @@ -70,7 +72,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) return drho; } -double Charge_Mixing::get_dkin(Charge* chr, const double nelec) +double Charge_Mixing::get_dkin(Charge* chr, const double nelec , const double omega) { if (!(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)) { @@ -99,9 +101,9 @@ double Charge_Mixing::get_dkin(Charge* chr, const double nelec) Parallel_Reduce::reduce_pool(dkin); #endif assert(nelec != 0); - assert(GlobalC::ucell.omega > 0); + assert(omega > 0); assert(this->rhopw->nxyz > 0); - dkin *= GlobalC::ucell.omega / static_cast(this->rhopw->nxyz); + dkin *= omega / static_cast(this->rhopw->nxyz); dkin /= nelec; ModuleBase::timer::tick("Charge_Mixing", "get_dkin"); @@ -121,7 +123,7 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c rhog2[is] = rho2 + is * this->rhopw->npw; } - static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2; + static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / tpiba2; static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI); double sum = 0.0; @@ -245,7 +247,7 @@ double Charge_Mixing::inner_product_recip_rho(std::complex* rho1, std::c #endif ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_rho"); - sum *= GlobalC::ucell.omega * 0.5; + sum *= omega * 0.5; delete[] rhog1; delete[] rhog2; @@ -285,7 +287,7 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_hartree"); ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree"); - static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2; + static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / tpiba2; static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI); double sum = 0.0; @@ -446,7 +448,7 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex* rhog1, s ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree"); - sum *= GlobalC::ucell.omega * 0.5; + sum *= omega * 0.5; return sum; } diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index 7d92e3959a..23dad6bd87 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -156,7 +156,7 @@ void H_TDDFT_pw::cal_v_space_length(std::vector& vext_space, int direc) ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length"); ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length"); - prepare(GlobalC::ucell, direc); + prepare(this->ucell_->G, direc); for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) { @@ -436,28 +436,26 @@ double H_TDDFT_pw::cal_v_time_heaviside() return vext_time; } -void H_TDDFT_pw::prepare(const UnitCell& cell, int& dir) +void H_TDDFT_pw::prepare(const ModuleBase::Matrix3& G, int& dir) { - if (dir == 1) - { - bvec[0] = cell.G.e11; - bvec[1] = cell.G.e12; - bvec[2] = cell.G.e13; - } - else if (dir == 2) - { - bvec[0] = cell.G.e21; - bvec[1] = cell.G.e22; - bvec[2] = cell.G.e23; - } - else if (dir == 3) - { - bvec[0] = cell.G.e31; - bvec[1] = cell.G.e32; - bvec[2] = cell.G.e33; - } - else + switch (dir) { + case 1: + bvec[0] = G.e11; + bvec[1] = G.e12; + bvec[2] = G.e13; + break; + case 2: + bvec[0] = G.e21; + bvec[1] = G.e22; + bvec[2] = G.e23; + break; + case 3: + bvec[0] = G.e31; + bvec[1] = G.e32; + bvec[2] = G.e33; + break; + default: ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!"); } bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index 4925da63bd..252b428fdc 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/source/module_elecstate/potentials/H_TDDFT_pw.h @@ -127,7 +127,7 @@ class H_TDDFT_pw : public PotBase //get ncut number for At integral static int check_ncut(int t_type); - void prepare(const UnitCell& cell, int& dir); + void prepare(const ModuleBase::Matrix3& G, int& dir); }; } // namespace elecstate diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp index b86f5dd64c..5a8ad4e6ce 100644 --- a/source/module_elecstate/test/charge_mixing_test.cpp +++ b/source/module_elecstate/test/charge_mixing_test.cpp @@ -48,10 +48,6 @@ InfoNonlocal::~InfoNonlocal() } #endif // mock class cell -namespace GlobalC -{ -UnitCell ucell; -}; /************************************************ * unit test of charge_mixing.cpp ***********************************************/ @@ -89,6 +85,7 @@ UnitCell ucell; class ChargeMixingTest : public ::testing::Test { public: + UnitCell ucell; ChargeMixingTest() { // Init pw_basis @@ -111,6 +108,8 @@ class ChargeMixingTest : public ::testing::Test PARAM.input.mixing_gg0_min = 0.1; PARAM.input.mixing_angle = -10.0; PARAM.input.mixing_dmr = false; + ucell.tpiba = 1.0; + ucell.omega = 2.0; } ModulePW::PW_Basis pw_basis; ModulePW::PW_Basis_Sup pw_dbasis; @@ -138,7 +137,9 @@ TEST_F(ChargeMixingTest, SetMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); EXPECT_EQ(CMtest.get_mixing_mode(), "broyden"); EXPECT_EQ(CMtest.get_mixing_beta(), 1.0); EXPECT_EQ(CMtest.get_mixing_ndim(), 1); @@ -161,7 +162,9 @@ TEST_F(ChargeMixingTest, SetMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); EXPECT_EQ(CMtest.mixing_mode, "plain"); EXPECT_EQ(CMtest.mixing_tau, true); @@ -177,7 +180,9 @@ TEST_F(ChargeMixingTest, SetMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), ""); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba);, ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("You'd better set mixing_beta to [0.0, 1.0]!")); @@ -194,7 +199,9 @@ TEST_F(ChargeMixingTest, SetMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), ""); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba);, ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("You'd better set mixing_beta_mag >= 0.0!")); @@ -212,7 +219,9 @@ TEST_F(ChargeMixingTest, SetMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), ""); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba);, ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon.")); } @@ -236,7 +245,9 @@ TEST_F(ChargeMixingTest, InitMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); PARAM.input.scf_thr_type= 1; CMtest.init_mixing(); @@ -261,7 +272,9 @@ TEST_F(ChargeMixingTest, InitMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); FUNC_TYPE = 3; CMtest.init_mixing(); EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.nrxx); @@ -277,7 +290,9 @@ TEST_F(ChargeMixingTest, InitMixingTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest.init_mixing(); EXPECT_EQ(CMtest.rho_mdata.length, 2 * pw_basis.nrxx); } @@ -295,7 +310,9 @@ TEST_F(ChargeMixingTest, InnerDotRealTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest.set_rhopw(&pw_basis, &pw_basis); PARAM.input.nspin = 4; @@ -321,7 +338,9 @@ TEST_F(ChargeMixingTest, InnerDotRealTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); PARAM.input.nspin = 4; // a simple sum for inner product @@ -349,7 +368,9 @@ TEST_F(ChargeMixingTest, InnerDotRecipSimpleTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest.set_rhopw(&pw_basis, &pw_basis); PARAM.input.nspin = 2; @@ -384,9 +405,18 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest) EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8); // RECIPROCAL NSPIN=1 - GlobalC::ucell.tpiba2 = 1.0; - GlobalC::ucell.omega = 2.0; - + CMtest.set_mixing(PARAM.input.mixing_mode, + PARAM.input.mixing_beta, + PARAM.input.mixing_ndim, + PARAM.input.mixing_gg0, + PARAM.input.mixing_tau, + PARAM.input.mixing_beta_mag, + PARAM.input.mixing_gg0_mag, + PARAM.input.mixing_gg0_min, + PARAM.input.mixing_angle, + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); PARAM.input.nspin = 1; std::vector> drhog1(pw_basis.npw); std::vector> drhog2(pw_basis.npw); @@ -459,7 +489,9 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); drhog1.resize(pw_basis.npw * 2); drhog2.resize(pw_basis.npw * 2); for (int i = 0; i < pw_basis.npw * 2; ++i) @@ -481,6 +513,18 @@ TEST_F(ChargeMixingTest, InnerDotRecipRhoTest) Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); PARAM.input.nspin = 1; + CMtest.set_mixing(PARAM.input.mixing_mode, + PARAM.input.mixing_beta, + PARAM.input.mixing_ndim, + PARAM.input.mixing_gg0, + PARAM.input.mixing_tau, + PARAM.input.mixing_beta_mag, + PARAM.input.mixing_gg0_mag, + PARAM.input.mixing_gg0_min, + PARAM.input.mixing_angle, + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); std::vector drhor1(pw_basis.nrxx); std::vector drhor2(pw_basis.nrxx); for (int i = 0; i < pw_basis.nrxx; ++i) @@ -492,8 +536,8 @@ TEST_F(ChargeMixingTest, InnerDotRecipRhoTest) EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8); // RECIPROCAL - GlobalC::ucell.tpiba2 = 1.0; - GlobalC::ucell.omega = 2.0; + ucell.tpiba = 1.0; + ucell.omega = 2.0; PARAM.input.nspin = 1; std::vector> drhog1(pw_basis.npw); @@ -548,7 +592,19 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest) { Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalC::ucell.tpiba = 1.0; + ucell.tpiba = 1.0; + CMtest.set_mixing(PARAM.input.mixing_mode, + PARAM.input.mixing_beta, + PARAM.input.mixing_ndim, + PARAM.input.mixing_gg0, + PARAM.input.mixing_tau, + PARAM.input.mixing_beta_mag, + PARAM.input.mixing_gg0_mag, + PARAM.input.mixing_gg0_min, + PARAM.input.mixing_angle, + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); // nspin = 1 PARAM.input.nspin = 1; std::complex* drhog = new std::complex[PARAM.input.nspin*pw_basis.npw]; @@ -678,8 +734,19 @@ TEST_F(ChargeMixingTest, KerkerScreenRealTest) { Charge_Mixing CMtest; CMtest.set_rhopw(&pw_basis, &pw_basis); - GlobalC::ucell.tpiba = 1.0; - + ucell.tpiba = 1.0; + CMtest.set_mixing(PARAM.input.mixing_mode, + PARAM.input.mixing_beta, + PARAM.input.mixing_ndim, + PARAM.input.mixing_gg0, + PARAM.input.mixing_tau, + PARAM.input.mixing_beta_mag, + PARAM.input.mixing_gg0_mag, + PARAM.input.mixing_gg0_min, + PARAM.input.mixing_angle, + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); // nspin = 1 PARAM.input.nspin = 1; double* drhor = new double[PARAM.input.nspin*pw_basis.nrxx]; @@ -815,7 +882,9 @@ TEST_F(ChargeMixingTest, MixRhoTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest_recip.init_mixing(); for(int i = 0 ; i < nspin * npw; ++i) { @@ -854,7 +923,9 @@ TEST_F(ChargeMixingTest, MixRhoTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest_real.init_mixing(); for(int i = 0 ; i < nspin * nrxx; ++i) { @@ -950,7 +1021,9 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest) PARAM.input.mixing_gg0_mag, PARAM.input.mixing_gg0_min, PARAM.input.mixing_angle, - PARAM.input.mixing_dmr); + PARAM.input.mixing_dmr, + ucell.omega, + ucell.tpiba); CMtest_recip.init_mixing(); for (int i = 0; i < nspin * npw; ++i) diff --git a/source/module_elecstate/test/elecstate_base_test.cpp b/source/module_elecstate/test/elecstate_base_test.cpp index 050a50b74a..cf9bb934ee 100644 --- a/source/module_elecstate/test/elecstate_base_test.cpp +++ b/source/module_elecstate/test/elecstate_base_test.cpp @@ -60,13 +60,14 @@ void ModulePW::PW_Basis::initgrids(double, ModuleBase::Matrix3, int, int, int) void ModulePW::PW_Basis::distribute_r() { } -void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*) +void Charge::set_rho_core(const UnitCell& , ModuleBase::ComplexMatrix const&, const bool*) { } void Charge::set_rho_core_paw() { } -void Charge::init_rho(elecstate::efermi&, +void Charge::init_rho(const UnitCell&, + elecstate::efermi&, ModuleBase::ComplexMatrix const&, ModuleSymmetry::Symmetry& symm, const void*, @@ -82,7 +83,24 @@ void Charge::renormalize_rho() void Charge::check_rho() { } - +UnitCell::UnitCell() +{ +} +UnitCell::~UnitCell() +{ +} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} /************************************************ * unit test of elecstate.cpp ***********************************************/ @@ -140,6 +158,7 @@ class ElecStateTest : public ::testing::Test protected: elecstate::MockElecState* elecstate; std::string output; + UnitCell ucell; void SetUp() { elecstate = new elecstate::MockElecState; @@ -249,7 +268,7 @@ TEST_F(ElecStateTest, InitSCF) ModuleBase::ComplexMatrix strucfac; elecstate->eferm = efermi; ModuleSymmetry::Symmetry symm; - EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, nullptr, symm)); + EXPECT_NO_THROW(elecstate->init_scf(istep, ucell,strucfac, nullptr, symm)); // delete elecstate->pot is done in the destructor of elecstate delete charge; } diff --git a/source/module_elecstate/test/elecstate_energy_test.cpp b/source/module_elecstate/test/elecstate_energy_test.cpp index 3d145ea187..28f01688a7 100644 --- a/source/module_elecstate/test/elecstate_energy_test.cpp +++ b/source/module_elecstate/test/elecstate_energy_test.cpp @@ -56,8 +56,24 @@ K_Vectors::K_Vectors() K_Vectors::~K_Vectors() { } - - +UnitCell::UnitCell() +{ +} +UnitCell::~UnitCell() +{ +} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} /*************************************************************** * unit test of functions in elecstate_energy.cpp ****************************************************************/ @@ -103,6 +119,7 @@ void ElecState::calculate_weights() class ElecStateEnergyTest : public ::testing::Test { protected: + UnitCell ucell; elecstate::MockElecState* elecstate; void SetUp() override { @@ -208,7 +225,7 @@ TEST_F(ElecStateEnergyTest, CalBandgap) TEST_F(ElecStateEnergyTest, CalBandgapUpDwTrivial) { - elecstate->cal_bandgap_updw(); + elecstate->cal_bandgap_updw(ucell); EXPECT_DOUBLE_EQ(elecstate->bandgap_up, 0.0); EXPECT_DOUBLE_EQ(elecstate->bandgap_dw, 0.0); } @@ -247,7 +264,7 @@ TEST_F(ElecStateEnergyTest, CalBandgapUpDw) } elecstate->eferm.ef_up = 0.5; elecstate->eferm.ef_dw = 2.1; - elecstate->cal_bandgap_updw(); + elecstate->cal_bandgap_updw(ucell); EXPECT_DOUBLE_EQ(elecstate->bandgap_up, 1.0); EXPECT_DOUBLE_EQ(elecstate->bandgap_dw, 0.5); } diff --git a/source/module_elecstate/test/elecstate_magnetism_test.cpp b/source/module_elecstate/test/elecstate_magnetism_test.cpp index b9266c7244..e826d533a6 100644 --- a/source/module_elecstate/test/elecstate_magnetism_test.cpp +++ b/source/module_elecstate/test/elecstate_magnetism_test.cpp @@ -31,10 +31,6 @@ Charge::~Charge() { } -double elecstate::get_ucell_omega() -{ - return 500.0; -} class MagnetismTest : public ::testing::Test { @@ -57,10 +53,6 @@ TEST_F(MagnetismTest, Magnetism) EXPECT_EQ(nullptr, magnetism->start_magnetization); } -TEST_F(MagnetismTest, GlobalInfo) -{ - EXPECT_EQ(500.0, elecstate::get_ucell_omega()); -} TEST_F(MagnetismTest, JudgeParallel) { @@ -91,7 +83,7 @@ TEST_F(MagnetismTest, ComputeMagnetizationS2) chr->rho[1][ir] = 1.01; } double* nelec_spin = new double[2]; - magnetism->compute_magnetization(chr->nrxx, chr->nxyz, chr->rho, nelec_spin); + magnetism->compute_magnetization(500.0, chr->nrxx, chr->nxyz, chr->rho, nelec_spin); EXPECT_DOUBLE_EQ(-0.5, magnetism->tot_magnetization); EXPECT_DOUBLE_EQ(0.5, magnetism->abs_magnetization); EXPECT_DOUBLE_EQ(4.75, nelec_spin[0]); @@ -125,7 +117,7 @@ TEST_F(MagnetismTest, ComputeMagnetizationS4) chr->rho[3][ir] = 1.00; } double* nelec_spin = new double[4]; - magnetism->compute_magnetization(chr->nrxx, chr->nxyz, chr->rho, nelec_spin); + magnetism->compute_magnetization(500.0, chr->nrxx, chr->nxyz, chr->rho, nelec_spin); EXPECT_DOUBLE_EQ(100.0, magnetism->abs_magnetization); EXPECT_DOUBLE_EQ(50.0*std::sqrt(2.0), magnetism->tot_magnetization_nc[0]); EXPECT_DOUBLE_EQ(50.0, magnetism->tot_magnetization_nc[1]); diff --git a/source/module_elecstate/test/elecstate_print_test.cpp b/source/module_elecstate/test/elecstate_print_test.cpp index 657ac999b8..859218386c 100644 --- a/source/module_elecstate/test/elecstate_print_test.cpp +++ b/source/module_elecstate/test/elecstate_print_test.cpp @@ -17,7 +17,12 @@ K_Vectors::K_Vectors() K_Vectors::~K_Vectors() { } - +UnitCell::UnitCell(){} +UnitCell::~UnitCell(){} +Magnetism::Magnetism(){} +Magnetism::~Magnetism(){} +InfoNonlocal::InfoNonlocal(){} +InfoNonlocal::~InfoNonlocal(){} /*************************************************************** * mock functions ****************************************************************/ @@ -33,11 +38,6 @@ void ElecState::calculate_weights() } // just for mock double Efield::etotefield = 1.1; double elecstate::Gatefield::etotgatefield = 2.2; -std::string tmp_vdw_method = "d2"; -std::string get_input_vdw_method() -{ - return tmp_vdw_method; -} double get_ucell_tot_magnetization() { return 1.1; @@ -58,11 +58,7 @@ double get_ucell_tot_magnetization_nc_z() { return 5.5; } -std::string tmp_ks_solver = "dav"; -std::string get_ks_solver_type() -{ - return tmp_ks_solver; -} + } // namespace elecstate Charge::Charge() @@ -94,6 +90,7 @@ class ElecStatePrintTest : public ::testing::Test std::ifstream ifs; std::ofstream ofs; K_Vectors* p_klist; + UnitCell ucell; void SetUp() { p_klist = new K_Vectors; @@ -117,6 +114,11 @@ class ElecStatePrintTest : public ::testing::Test elecstate.wg(0, 1) = 0.2; elecstate.wg(1, 0) = 0.3; elecstate.wg(1, 1) = 0.4; + ucell.magnet.tot_magnetization = 1.1; + ucell.magnet.abs_magnetization = 2.2; + ucell.magnet.tot_magnetization_nc[0] = 3.3; + ucell.magnet.tot_magnetization_nc[1] = 4.4; + ucell.magnet.tot_magnetization_nc[2] = 5.5; } void TearDown() { @@ -251,38 +253,38 @@ TEST_F(ElecStatePrintTest, PrintEtot) std::vector vdw_methods = {"d2", "d3_0", "d3_bj"}; for (int i = 0; i < vdw_methods.size(); i++) { - elecstate::tmp_vdw_method = vdw_methods[i]; - elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, false); + PARAM.input.vdw_method = vdw_methods[i]; + elecstate.print_etot(ucell,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, false); } // iteration of different ks_solver std::vector ks_solvers = {"cg", "lapack", "genelpa", "dav", "scalapack_gvx", "cusolver"}; for (int i = 0; i < ks_solvers.size(); i++) { - elecstate::tmp_ks_solver = ks_solvers[i]; + PARAM.input.ks_solver = ks_solvers[i]; testing::internal::CaptureStdout(); - elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); + elecstate.print_etot(ucell,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); output = testing::internal::GetCapturedStdout(); - if (elecstate::tmp_ks_solver == "cg") + if (PARAM.input.ks_solver == "cg") { EXPECT_THAT(output, testing::HasSubstr("CG")); } - else if (elecstate::tmp_ks_solver == "lapack") + else if (PARAM.input.ks_solver == "lapack") { EXPECT_THAT(output, testing::HasSubstr("LA")); } - else if (elecstate::tmp_ks_solver == "genelpa") + else if (PARAM.input.ks_solver == "genelpa") { EXPECT_THAT(output, testing::HasSubstr("GE")); } - else if (elecstate::tmp_ks_solver == "dav") + else if (PARAM.input.ks_solver == "dav") { EXPECT_THAT(output, testing::HasSubstr("DA")); } - else if (elecstate::tmp_ks_solver == "scalapack_gvx") + else if (PARAM.input.ks_solver == "scalapack_gvx") { EXPECT_THAT(output, testing::HasSubstr("GV")); } - else if (elecstate::tmp_ks_solver == "cusolver") + else if (PARAM.input.ks_solver == "cusolver") { EXPECT_THAT(output, testing::HasSubstr("CU")); } @@ -294,7 +296,7 @@ TEST_F(ElecStatePrintTest, PrintEtot) EXPECT_THAT(str, testing::HasSubstr("Error Threshold = 0.1")); EXPECT_THAT(str, testing::HasSubstr("E_KohnSham")); EXPECT_THAT(str, testing::HasSubstr("E_vdwD2")); - EXPECT_THAT(str, testing::HasSubstr("E_vdwD3")); + // EXPECT_THAT(str, testing::HasSubstr("E_vdwD3")); EXPECT_THAT(str, testing::HasSubstr("E_sol_el")); EXPECT_THAT(str, testing::HasSubstr("E_sol_cav")); EXPECT_THAT(str, testing::HasSubstr("E_efield")); @@ -328,7 +330,7 @@ TEST_F(ElecStatePrintTest, PrintEtot2) PARAM.input.basis_type = "pw"; PARAM.input.scf_nmax = 100; - elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); + elecstate.print_etot(ucell,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); GlobalV::ofs_running.close(); ifs.open("test.dat", std::ios::in); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); @@ -364,7 +366,7 @@ TEST_F(ElecStatePrintTest, PrintEtotColorS2) PARAM.input.out_bandgap = true; PARAM.input.nspin = 2; GlobalV::MY_RANK = 0; - elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); + elecstate.print_etot(ucell,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); } TEST_F(ElecStatePrintTest, PrintEtotColorS4) @@ -389,5 +391,5 @@ TEST_F(ElecStatePrintTest, PrintEtotColorS4) PARAM.input.nspin = 4; 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); + elecstate.print_etot(ucell,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print); } \ No newline at end of file diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp index 85ba4d1b24..9d85695328 100644 --- a/source/module_elecstate/test/elecstate_pw_test.cpp +++ b/source/module_elecstate/test/elecstate_pw_test.cpp @@ -141,13 +141,14 @@ K_Vectors::~K_Vectors() { } -void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*) +void Charge::set_rho_core(const UnitCell& ,ModuleBase::ComplexMatrix const&, const bool*) { } void Charge::set_rho_core_paw() { } -void Charge::init_rho(elecstate::efermi&, +void Charge::init_rho(const UnitCell&, + elecstate::efermi&, ModuleBase::ComplexMatrix const&, ModuleSymmetry::Symmetry& symm, const void*, diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index a79ab77493..582e63b137 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -108,7 +108,9 @@ void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para PARAM.inp.mixing_gg0_mag, PARAM.inp.mixing_gg0_min, PARAM.inp.mixing_angle, - PARAM.inp.mixing_dmr); + PARAM.inp.mixing_dmr, + ucell.omega, + ucell.tpiba); p_chgmix->init_mixing(); /// PAW Section @@ -368,7 +370,7 @@ void ESolver_KS::hamilt2density(UnitCell& ucell, const int istep, con // double drho = this->estate.caldr2(); // EState should be used after it is constructed. - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec,ucell.omega); hsolver_error = 0.0; if (iter == 1 && PARAM.inp.calculation != "nscf") { @@ -390,7 +392,7 @@ void ESolver_KS::hamilt2density(UnitCell& ucell, const int istep, con this->hamilt2density_single(ucell, istep, iter, diag_ethr); - drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec); + drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec ,ucell.omega); hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type, PARAM.inp.esolver_type, @@ -540,7 +542,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i } else { - this->pelec->cal_bandgap_updw(); + this->pelec->cal_bandgap_updw(ucell); } } @@ -550,7 +552,8 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i } // compute magnetization, only for LSDA(spin==2) - ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, + ucell.magnet.compute_magnetization(ucell.omega, + this->pelec->charge->nrxx, this->pelec->charge->nxyz, this->pelec->charge->rho, this->pelec->nelec_spin.data()); @@ -673,9 +676,9 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i double dkin = 0.0; // for meta-GGA if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec); + dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec,ucell.omega); } - this->pelec->print_etot(this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr); + this->pelec->print_etot(ucell,this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr); // Json, need to be moved to somewhere else #ifdef __RAPIDJSON diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index e6a4ab6632..3e48e4aff0 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -558,7 +558,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const } // mohan update 2012-06-05 - this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell); // mohan move it outside 2011-01-13 // first need to calculate the weight according to @@ -746,7 +746,7 @@ void ESolver_KS_LCAO::hamilt2density_single(UnitCell& ucell, int istep, } // 12) calculate delta energy - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); } //------------------------------------------------------------------------------ @@ -1052,7 +1052,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) occ_number_ks(ik, inb) /= this->kv.wk[ik]; } } - this->rdmft_solver.update_elec(occ_number_ks, *(this->psi)); + this->rdmft_solver.update_elec(ucell,occ_number_ks, *(this->psi)); //! initialize the gradients of Etotal with respect to occupation numbers and wfc, //! and set all elements to 0. diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 880daa8512..d29bb6a6fc 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -139,7 +139,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(UnitCell& ucell, const int ist } // (7) calculate delta energy - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); } void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep, int& iter) diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index 1060a80d29..e79346cd6d 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -166,7 +166,7 @@ namespace ModuleESolver // deband is calculated from "output" charge density calculated // in sum_band // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))' - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density_single"); } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 33be890791..34532835b7 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -302,7 +302,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) elecstate::cal_ux(ucell); //! calculate the total local pseudopotential in real space - this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc); + this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) @@ -403,7 +403,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const } // mohan move harris functional to here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) - this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell); } // Temporary, it should be replaced by hsolver later. @@ -462,7 +462,7 @@ void ESolver_KS_PW::hamilt2density_single(UnitCell& ucell, // deband is calculated from "output" charge density calculated // in sum_band // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))' - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density_single"); } diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 32d4a90b2e..c3f8e794c2 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -120,7 +120,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) this->init_elecstate(ucell); // calculate the total local pseudopotential in real space - this->pelec->init_scf(0, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs + this->pelec->init_scf(0, ucell,sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -255,7 +255,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) GlobalV::ofs_warning); } - this->pelec->init_scf(istep, sf.strucFac, locpp.numeric, ucell.symm); + this->pelec->init_scf(istep, ucell,sf.strucFac, locpp.numeric, ucell.symm); // calculate ewald energy this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac); diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index fc18e8e302..8855592d30 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -188,7 +188,7 @@ void ESolver_SDFT_PW::hamilt2density_single(UnitCell& ucell, int iste { srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm); } - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); } else { diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index b02c92729d..0dcc7744db 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -264,7 +264,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) } #endif // __EXX - this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm); + this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm); //! output the initial charge density if (PARAM.inp.out_chg[0] == 2) diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 50012599c1..bdea01d17f 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -254,7 +254,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) elecstate::cal_ux(ucell); // pelec should be initialized before these calculations - this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm); + this->pelec->init_scf(istep,ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm); // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 389119a17a..fbfdb29fd8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -35,7 +35,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const bool isstress, const bool istestf, const bool istests, - const UnitCell& ucell, + UnitCell& ucell, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -851,7 +851,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, @@ -878,7 +878,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. //-------------------------------------------------------- @@ -992,7 +992,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, @@ -1024,7 +1024,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(ucell,sigmacc,rhopw, &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 927d08790e..28b52a9975 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -33,7 +33,7 @@ class Force_Stress_LCAO const bool isstress, const bool istestf, const bool istests, - const UnitCell& ucell, + UnitCell& ucell, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -64,7 +64,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, @@ -102,7 +102,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/hamilt_lcaodft/LCAO_hamilt.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp index 42317c2df7..16221587cf 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp @@ -22,6 +22,7 @@ template void sparse_format::cal_HR_exx( + const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, @@ -36,15 +37,15 @@ void sparse_format::cal_HR_exx( const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha; std::map> atoms_pos; - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) { + for (int iat = 0; iat < ucell.nat; ++iat) { atoms_pos[iat] = RI_Util::Vector3_to_array3( - GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]] - .tau[GlobalC::ucell.iat2ia[iat]]); + ucell.atoms[ucell.iat2it[iat]] + .tau[ucell.iat2ia[iat]]); } const std::array, 3> latvec - = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), // too bad to use GlobalC here, - RI_Util::Vector3_to_array3(GlobalC::ucell.a2), - RI_Util::Vector3_to_array3(GlobalC::ucell.a3)}; + = {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here, + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3)}; const std::array Rs_period = {nmp[0], nmp[1], nmp[2]}; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index f708ad1310..bba142d8f4 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -223,6 +223,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, dftu = new OperatorDFTU>(this->hsk, kv->kvec_d, this->hR, // no explicit call yet + ucell, this->kv->isk); } else @@ -365,6 +366,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, dftu = new OperatorDFTU>(this->hsk, kv->kvec_d, this->hR, // no explicit call yet + ucell, this->kv->isk); } else diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.cpp index 924dcf1764..c4eade40f9 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.cpp @@ -27,7 +27,7 @@ void OperatorDFTU>::contributeHk(int ik) ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_real(ik, &eff_pot[0], isk, this->hsk->get_sk()); + GlobalC::dftu.cal_eff_pot_mat_real(ik, ucell, &eff_pot[0], isk, this->hsk->get_sk()); double* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) @@ -45,7 +45,7 @@ void OperatorDFTU, double>>::contributeHk(int ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + GlobalC::dftu.cal_eff_pot_mat_complex(ik, ucell,&eff_pot[0], isk, this->hsk->get_sk()); std::complex* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) @@ -63,7 +63,7 @@ void OperatorDFTU, std::complex>>::con ModuleBase::timer::tick("OperatorDFTU", "contributeHk"); // Effective potential of DFT+U is added to total Hamiltonian here; Quxin adds on 20201029 std::vector> eff_pot(this->hsk->get_pv()->nloc); - GlobalC::dftu.cal_eff_pot_mat_complex(ik, &eff_pot[0], isk, this->hsk->get_sk()); + GlobalC::dftu.cal_eff_pot_mat_complex(ik, ucell,&eff_pot[0], isk, this->hsk->get_sk()); std::complex* hk = this->hsk->get_hk(); for (int irc = 0; irc < this->hsk->get_pv()->nloc; irc++) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h index de40bce4b3..5ba2745597 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h @@ -23,8 +23,9 @@ class OperatorDFTU> : public OperatorLCAO OperatorDFTU>(HS_Matrix_K* hsk_in, const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, + const UnitCell& ucell_in, const std::vector& isk_in) - : isk(isk_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) + : ucell(ucell_in),isk(isk_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_dftu; } @@ -35,6 +36,8 @@ class OperatorDFTU> : public OperatorLCAO private: + const UnitCell& ucell; + bool HR_fixed_done = false; const std::vector& isk; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h index a9b7b7156f..775c60fdaf 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h @@ -9,6 +9,7 @@ #include #include +#include "module_cell/unitcell.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp" #include "module_basis/module_ao/parallel_orbitals.h" @@ -17,6 +18,7 @@ namespace sparse_format template void cal_HR_exx( + const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index d845c3c8f2..c2febdbb25 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -7,7 +7,8 @@ #include "spar_exx.h" #include "spar_u.h" -void sparse_format::cal_HSR(const Parallel_Orbitals& pv, +void sparse_format::cal_HSR(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const int& current_spin, @@ -76,6 +77,7 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, if (PARAM.inp.dft_plus_u == 2) { if (nspin == 1 || nspin == 2) { cal_HR_dftu(pv, + ucell, HS_Arrays.all_R_coor, HS_Arrays.SR_sparse, HS_Arrays.HR_sparse, @@ -83,6 +85,7 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, sparse_thr); } else if (nspin == 4) { cal_HR_dftu_soc(pv, + ucell, HS_Arrays.all_R_coor, HS_Arrays.SR_soc_sparse, HS_Arrays.HR_soc_sparse, @@ -99,7 +102,9 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, if (GlobalC::exx_info.info_global.cal_exx) { if (Hexxd && GlobalC::exx_info.info_ri.real_number) { - sparse_format::cal_HR_exx(pv, + sparse_format::cal_HR_exx( + ucell, + pv, HS_Arrays, current_spin, sparse_thr, @@ -108,7 +113,9 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, } else if (Hexxc && !GlobalC::exx_info.info_ri.real_number) { - sparse_format::cal_HR_exx(pv, + sparse_format::cal_HR_exx( + ucell, + pv, HS_Arrays, current_spin, sparse_thr, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index 497586574e..a1118f994e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -5,7 +5,8 @@ namespace sparse_format { using TAC = std::pair>; - void cal_HSR(const Parallel_Orbitals& pv, + void cal_HSR(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const int& current_spin, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.cpp index 4515ea185c..031e1e671e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.cpp @@ -7,6 +7,7 @@ void sparse_format::cal_HR_dftu( const Parallel_Orbitals &pv, + const UnitCell& ucell, std::set> &all_R_coor, std::map, std::map>> &SR_sparse, std::map, std::map>> *HR_sparse, @@ -74,7 +75,7 @@ void sparse_format::cal_HR_dftu( } } - GlobalC::dftu.cal_eff_pot_mat_R_double(current_spin, SR_tmp, HR_tmp); + GlobalC::dftu.cal_eff_pot_mat_R_double(current_spin, ucell,SR_tmp, HR_tmp); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { @@ -129,6 +130,7 @@ void sparse_format::cal_HR_dftu( void sparse_format::cal_HR_dftu_soc( const Parallel_Orbitals &pv, + const UnitCell &ucell, std::set> &all_R_coor, std::map, std::map>>> &SR_soc_sparse, std::map, std::map>>> &HR_soc_sparse, @@ -194,7 +196,7 @@ void sparse_format::cal_HR_dftu_soc( } } - GlobalC::dftu.cal_eff_pot_mat_R_complex_double(current_spin, SR_soc_tmp, HR_soc_tmp); + GlobalC::dftu.cal_eff_pot_mat_R_complex_double(current_spin, ucell, SR_soc_tmp, HR_soc_tmp); for (int i = 0; i < PARAM.globalv.nlocal; ++i) { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.h index e26d808b5a..bdf48ef15d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_u.h @@ -12,6 +12,7 @@ namespace sparse_format void cal_HR_dftu( const Parallel_Orbitals &pv, + const UnitCell &ucell, std::set> &all_R_coor, std::map, std::map>> &SR_sparse, std::map, std::map>> *HR_sparse, @@ -20,6 +21,7 @@ namespace sparse_format void cal_HR_dftu_soc( const Parallel_Orbitals &pv, + const UnitCell &ucell, std::set> &all_R_coor, std::map, std::map>>> &SR_soc_sparse, std::map, std::map>>> &HR_soc_sparse, diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index 9543ae6e55..da11baeaac 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -77,10 +77,10 @@ class DFTU // For calculating contribution to Hamiltonian matrices //============================================================= public: - void cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::vector& isk, const std::complex* sk); - void cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk); - void cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR); - void cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR); + void cal_eff_pot_mat_complex(const int ik, const UnitCell& ucell, std::complex* eff_pot, const std::vector& isk, const std::complex* sk); + void cal_eff_pot_mat_real(const int ik, const UnitCell& ucell, double* eff_pot, const std::vector& isk, const double* sk); + void cal_eff_pot_mat_R_double(const int ispin,const UnitCell& ucell, double* SR, double* HR); + void cal_eff_pot_mat_R_complex_double(const int ispin, const UnitCell& ucell, std::complex* SR, std::complex* HR); //============================================================= // In dftu_occup.cpp @@ -122,8 +122,14 @@ class DFTU // for both Hamiltonian and force/stress //============================================================= - void cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU); - void cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU); + void cal_VU_pot_mat_complex(const int spin, + const bool newlocale, + const UnitCell& ucell, + std::complex* VU); + void cal_VU_pot_mat_real(const int spin, + const bool newlocale, + const UnitCell& ucell, + double* VU); double get_onebody_eff_pot(const int T, const int iat, diff --git a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp index b24aa09865..923db1e7ac 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp @@ -110,7 +110,7 @@ void DFTU::force_stress(const UnitCell& ucell, double* VU = new double[pv.nloc]; - this->cal_VU_pot_mat_real(spin, false, VU); + this->cal_VU_pot_mat_real(spin, false, ucell, VU); const std::vector>& dmk = dynamic_cast*>(pelec)->get_DM()->get_DMK_vector(); @@ -174,7 +174,7 @@ void DFTU::force_stress(const UnitCell& ucell, std::complex* VU = new std::complex[pv.nloc]; - this->cal_VU_pot_mat_complex(spin, false, VU); + this->cal_VU_pot_mat_complex(spin, false, ucell, VU); const std::vector>>& dmk = dynamic_cast>*>(pelec) diff --git a/source/module_hamilt_lcao/module_dftu/dftu_hamilt.cpp b/source/module_hamilt_lcao/module_dftu/dftu_hamilt.cpp index c350aef529..19b05f90a9 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_hamilt.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_hamilt.cpp @@ -7,7 +7,11 @@ namespace ModuleDFTU { -void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::vector& isk, const std::complex* sk) +void DFTU::cal_eff_pot_mat_complex(const int ik, + const UnitCell& ucell, + std::complex* eff_pot, + const std::vector& isk, + const std::complex* sk) { ModuleBase::TITLE("DFTU", "cal_eff_pot_mat"); ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); @@ -31,7 +35,7 @@ void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, const std::complex zero = 0.0; std::vector> VU(this->paraV->nloc); - this->cal_VU_pot_mat_complex(spin, true, &VU[0]); + this->cal_VU_pot_mat_complex(spin, true, ucell,&VU[0]); #ifdef __MPI pzgemm_(&transN, &transN, @@ -58,7 +62,11 @@ void DFTU::cal_eff_pot_mat_complex(const int ik, std::complex* eff_pot, return; } -void DFTU::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector& isk, const double* sk) +void DFTU::cal_eff_pot_mat_real(const int ik, + const UnitCell& ucell, + double* eff_pot, + const std::vector& isk, + const double* sk) { ModuleBase::TITLE("DFTU", "cal_eff_pot_mat"); ModuleBase::timer::tick("DFTU", "cal_eff_pot_mat"); @@ -80,7 +88,7 @@ void DFTU::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector double alpha = 1.0, beta = 0.0, half = 0.5, one = 1.0; std::vector VU(this->paraV->nloc); - this->cal_VU_pot_mat_real(spin, 1, &VU[0]); + this->cal_VU_pot_mat_real(spin, 1, ucell, &VU[0]); #ifdef __MPI pdgemm_(&transN, &transN, @@ -107,14 +115,14 @@ void DFTU::cal_eff_pot_mat_real(const int ik, double* eff_pot, const std::vector return; } -void DFTU::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) +void DFTU::cal_eff_pot_mat_R_double(const int ispin, const UnitCell& ucell,double* SR, double* HR) { const char transN = 'N', transT = 'T'; const int one_int = 1; const double alpha = 1.0, beta = 0.0, one = 1.0, half = 0.5; std::vector VU(this->paraV->nloc); - this->cal_VU_pot_mat_real(ispin, 1, &VU[0]); + this->cal_VU_pot_mat_real(ispin, 1, ucell, &VU[0]); #ifdef __MPI pdgemm_(&transN, &transN, @@ -137,14 +145,14 @@ void DFTU::cal_eff_pot_mat_R_double(const int ispin, double* SR, double* HR) return; } -void DFTU::cal_eff_pot_mat_R_complex_double(const int ispin, std::complex* SR, std::complex* HR) +void DFTU::cal_eff_pot_mat_R_complex_double(const int ispin, const UnitCell& ucell, std::complex* SR, std::complex* HR) { const char transN = 'N', transT = 'T'; const int one_int = 1; const std::complex zero = 0.0, one = 1.0, half = 0.5; std::vector> VU(this->paraV->nloc); - this->cal_VU_pot_mat_complex(ispin, 1, &VU[0]); + this->cal_VU_pot_mat_complex(ispin, 1, ucell, &VU[0]); #ifdef __MPI pzgemm_(&transN, &transN, diff --git a/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp b/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp index 96582ee6aa..000642c09f 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp @@ -6,30 +6,33 @@ namespace ModuleDFTU { -void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::complex* VU) +void DFTU::cal_VU_pot_mat_complex(const int spin, + const bool newlocale, + const UnitCell& ucell, + std::complex* VU) { 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 < ucell.ntype; ++it) { if (PARAM.inp.orbital_corr[it] == -1) { continue; -} - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + } + for (int ia = 0; ia < 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 = ucell.itia2iat(it, ia); + for (int L = 0; L <= ucell.atoms[it].nwl; L++) { 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 < ucell.atoms[it].l_nchi[L]; n++) { if (n != 0) { continue; -} + } for (int m1 = 0; m1 < 2 * L + 1; m1++) { @@ -38,7 +41,7 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) { continue; -} + } for (int m2 = 0; m2 < 2 * L + 1; m2++) { @@ -48,7 +51,7 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); if (nu < 0) { continue; -} + } int m1_all = m1 + (2 * L + 1) * ipol1; int m2_all = m2 + (2 * L + 1) * ipol2; @@ -67,30 +70,33 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com return; } -void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) +void DFTU::cal_VU_pot_mat_real(const int spin, + const bool newlocale, + const UnitCell& ucell, + 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 < ucell.ntype; ++it) { if (PARAM.inp.orbital_corr[it] == -1) { continue; -} - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + } + for (int ia = 0; ia < 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 = ucell.itia2iat(it, ia); + for (int L = 0; L <= ucell.atoms[it].nwl; L++) { 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 < ucell.atoms[it].l_nchi[L]; n++) { if (n != 0) { continue; -} + } for (int m1 = 0; m1 < 2 * L + 1; m1++) { @@ -99,7 +105,7 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]); if (mu < 0) { continue; -} + } for (int m2 = 0; m2 < 2 * L + 1; m2++) { for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++) @@ -108,7 +114,7 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU) = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]); if (nu < 0) { continue; -} + } int m1_all = m1 + (2 * L + 1) * ipol1; int m2_all = m2 + (2 * L + 1) * ipol2; @@ -161,7 +167,7 @@ double DFTU::get_onebody_eff_pot(const int T, * (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 { @@ -169,7 +175,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 @@ -182,7 +188,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 { @@ -190,7 +196,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..cf22a3879b 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(ucell,sigmaxcc, this->rhopw, 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 f9c6a63556..940ed36e0c 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, @@ -162,7 +162,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 c23f24f53e..a1b5fe4288 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/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index 12162d88e4..f8c688de1c 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -33,7 +33,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, ModulePW::PW_Basis* rho_basis, const Charge* const chr, const bool* numeric, - UnitCell& ucell_in) + UnitCell& ucell_in) { ModuleBase::TITLE("Forces", "cal_force_cc"); // recalculate the exchange-correlation potential. diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h index 6d5ee5581e..c1399683ac 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h @@ -112,7 +112,8 @@ class Stress_Func ModulePW::PW_Basis* rho_basis); // used in local pseudopotential stress // 5) the stress from the non-linear core correction (if any) - void stress_cc(ModuleBase::matrix& sigma, + void stress_cc(UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf, const bool is_pw, 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..ab5ef87623 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -13,7 +13,8 @@ //NLCC term, need to be tested template -void Stress_Func::stress_cc(ModuleBase::matrix& sigma, +void Stress_Func::stress_cc(UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf, const bool is_pw, @@ -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,17 @@ 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.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); @@ -130,15 +131,15 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, local_sigmadiag = conj(psic[ig]) * p_sf->strucFac(nt, ig) * rhocg[rho_basis->ig2igg[ig]]; } else { local_sigmadiag = conj(psic[ig]) * p_sf->strucFac(nt, ig) * rhocg[rho_basis->ig2igg[ig]] * fact; -} + } sigmadiag += local_sigmadiag.real(); } 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.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); @@ -154,15 +155,15 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, for(int ig = 0;ig< rho_basis->npw;ig++) { const FPTYPE norm_g = sqrt(rho_basis->gg[ig]); - if(norm_g < 1e-4) { continue; -} + if(norm_g < 1e-4) {continue; + } for (int l = 0; l < 3; l++) { for (int m = 0;m< 3;m++) { 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 @@ -254,7 +255,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 / ucell->omega; igl0 = 1; } else @@ -273,7 +274,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] * ucell->tpiba2); } double *r_d = nullptr; @@ -298,12 +299,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,ucell->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,ucell->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 2a6925c912..6d95248c11 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp @@ -93,7 +93,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(ucell,sigmaxcc, rho_basis, 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..763b874947 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(ucell_in,sigmaxcc, rho_basis, 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_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h index 2ee381ad48..3ec0817825 100644 --- a/source/module_hsolver/test/hsolver_supplementary_mock.h +++ b/source/module_hsolver/test/hsolver_supplementary_mock.h @@ -43,6 +43,7 @@ void ElecState::print_eigenvalue(std::ofstream& ofs) } void ElecState::init_scf(const int istep, + const UnitCell& ucell, const ModuleBase::ComplexMatrix& strucfac, const bool*, ModuleSymmetry::Symmetry&, diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index 11cd2caf81..b182559e1e 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -45,7 +45,7 @@ void output_mat_sparse(const bool& out_mat_hsR, //! generate a file containing the Hamiltonian and S(overlap) matrices if (out_mat_hsR) { - output_HSR(istep, v_eff, pv, HS_Arrays, grid, kv, p_ham); + output_HSR(istep, ucell, v_eff, pv, HS_Arrays, grid, kv, p_ham); } //! generate a file containing the kinetic energy matrix diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 5b59f2ee24..549ad7cadb 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -13,6 +13,7 @@ // If the absolute value of the matrix element is less than or equal to the // 'sparse_thr', it will be ignored. void ModuleIO::output_HSR(const int& istep, + const UnitCell& ucell, const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, @@ -37,7 +38,7 @@ void ModuleIO::output_HSR(const int& istep, if (nspin == 1 || nspin == 4) { const int spin_now = 0; // jingan add 2021-6-4, modify 2021-12-2 - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif @@ -47,7 +48,7 @@ void ModuleIO::output_HSR(const int& istep, int spin_now = 1; // save HR of spin down first (the current spin always be down) - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif @@ -61,7 +62,7 @@ void ModuleIO::output_HSR(const int& istep, spin_now = 0; } - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h index 6c9bd46768..6473c369ec 100644 --- a/source/module_io/write_HS_R.h +++ b/source/module_io/write_HS_R.h @@ -12,6 +12,7 @@ namespace ModuleIO { using TAC = std::pair>; void output_HSR(const int& istep, + const UnitCell& ucell, const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp index d6fad87dad..e400e7ae8e 100644 --- a/source/module_io/write_vxc.hpp +++ b/source/module_io/write_vxc.hpp @@ -252,7 +252,7 @@ void write_Vxc(const int nspin, &vxcs_R_ao[0]/*for paraV*/, kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); std::vector> e_orb_exx; // orbital energy (EXX) #endif - hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, kv.isk); + hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, ucell,kv.isk); // 4. calculate and write the MO-matrix Exc Parallel_2D p2d; diff --git a/source/module_rdmft/rdmft.cpp b/source/module_rdmft/rdmft.cpp index e50d63fd45..ee1e0ee142 100644 --- a/source/module_rdmft/rdmft.cpp +++ b/source/module_rdmft/rdmft.cpp @@ -333,7 +333,7 @@ void RDMFT::cal_Energy(const int cal_type) } else { - this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + this->pelec->f_en.deband = this->pelec->cal_delta_eband(*ucell); E_descf = pelec->f_en.descf = 0.0; this->pelec->cal_energies(2); Etotal = this->pelec->f_en.etot; diff --git a/source/module_rdmft/rdmft.h b/source/module_rdmft/rdmft.h index 85e07c0613..0b381ce073 100644 --- a/source/module_rdmft/rdmft.h +++ b/source/module_rdmft/rdmft.h @@ -90,7 +90,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 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(); @@ -122,7 +122,7 @@ class RDMFT //! get the total Hamilton in k-space void cal_Hk_Hpsi(); - void update_charge(); + void update_charge(UnitCell& ucell_in); private: diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index a581910d6d..68c543b7e3 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -55,7 +55,7 @@ void RDMFT::cal_V_TV() hsk_TV, kv->kvec_d, HR_TV, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, two_center_bundle->kinetic_orb.get() @@ -65,7 +65,7 @@ void RDMFT::cal_V_TV() hsk_TV, kv->kvec_d, HR_TV, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, two_center_bundle->overlap_orb_beta.get() @@ -79,7 +79,7 @@ void RDMFT::cal_V_TV() kv->kvec_d, this->pelec->pot, HR_TV, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, @@ -98,7 +98,7 @@ void RDMFT::cal_V_TV() kv->kvec_d, this->pelec->pot, HR_TV, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, @@ -131,7 +131,7 @@ void RDMFT::cal_V_hartree() kv->kvec_d, this->pelec->pot, HR_hartree, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, @@ -151,7 +151,7 @@ void RDMFT::cal_V_hartree() kv->kvec_d, this->pelec->pot, HR_hartree, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, @@ -182,7 +182,7 @@ void RDMFT::cal_V_XC() // elecstate::DensityMatrix DM_test(ParaV, nspin, kv->kvec_d, nk_total); // elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test); - // DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + // DM_test.init_DMR(&GlobalC::GridD, ucell); // DM_test.cal_DMR(); // // compare DM_XC and DM get in update_charge(or ABACUS) @@ -215,7 +215,7 @@ void RDMFT::cal_V_XC() kv->kvec_d, this->pelec->pot, HR_dft_XC, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, @@ -237,7 +237,7 @@ void RDMFT::cal_V_XC() kv->kvec_d, this->pelec->pot, HR_dft_XC, - &GlobalC::ucell, + ucell, orb->cutoffs(), &GlobalC::GridD, nspin, diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index 3e6bbad183..020b648360 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -45,7 +45,10 @@ void RDMFT::update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis template -void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in) +void RDMFT::update_elec(UnitCell& ucell_in, + const ModuleBase::matrix& occ_number_in, + const psi::Psi& wfc_in, + const Charge* charge_in) { // update occ_number, wg, wk_fun_occNum occ_number = (occ_number_in); @@ -67,7 +70,7 @@ void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const p } // update charge - this->update_charge(); + this->update_charge(ucell_in); // "default" = "pbe" // if( !only_exx_type || this->cal_E_type != 1 ) @@ -87,14 +90,14 @@ void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const p // this code is copying from function ElecStateLCAO::psiToRho(), in elecstate_lcao.cpp template -void RDMFT::update_charge() +void RDMFT::update_charge(UnitCell& ucell_in) { 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(&GlobalC::GridD, &GlobalC::ucell); + DM_gamma_only.init_DMR(&GlobalC::GridD, ucell); DM_gamma_only.cal_DMR(); for (int is = 0; is < nspin; is++) @@ -124,7 +127,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(&GlobalC::GridD, &GlobalC::ucell); + DM.init_DMR(&GlobalC::GridD, ucell); DM.cal_DMR(); for (int is = 0; is < nspin; is++) @@ -156,7 +159,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_in.symm); } }