diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp index 48d25ed298..696a3a75a2 100644 --- a/source/module_elecstate/elecstate_pw.cpp +++ b/source/module_elecstate/elecstate_pw.cpp @@ -263,7 +263,7 @@ void ElecStatePW::add_usrho(const psi::Psi& psi) // get |beta> if (this->ppcell->nkb > 0) { - this->ppcell->getvnl(this->ctx, ik, this->vkb); + this->ppcell->getvnl(this->ctx, *ucell,ik, this->vkb); } // becp = diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp index c8692a01e0..85ba4d1b24 100644 --- a/source/module_elecstate/test/elecstate_pw_test.cpp +++ b/source/module_elecstate/test/elecstate_pw_test.cpp @@ -115,12 +115,14 @@ std::complex* pseudopot_cell_vnl::get_vkb_data() const } template <> void pseudopot_cell_vnl::getvnl(base_device::DEVICE_CPU*, + const UnitCell&, int const&, std::complex*) const { } template <> void pseudopot_cell_vnl::getvnl(base_device::DEVICE_CPU*, + const UnitCell&, int const&, std::complex*) const { diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 776e837d9a..10c3fc2b9b 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -281,7 +281,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rhod->collect_uniqgg(); } - this->p_locpp->init_vloc(this->pw_rhod); + this->p_locpp->init_vloc(ucell,this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); this->pelec->omega = ucell.omega; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index e2c83d8cf6..31d08fb71f 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -197,7 +197,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } // 8) initialize ppcell - this->ppcell.init_vloc(this->pw_rho); + this->ppcell.init_vloc(ucell,this->pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); // 9) inititlize the charge density diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp index e477ffecde..1060a80d29 100644 --- a/source/module_esolver/esolver_ks_lcaopw.cpp +++ b/source/module_esolver/esolver_ks_lcaopw.cpp @@ -62,9 +62,9 @@ namespace ModuleESolver } template - void ESolver_KS_LIP::allocate_hamilt() + void ESolver_KS_LIP::allocate_hamilt(const UnitCell& ucell) { - this->p_hamilt = new hamilt::HamiltLIP(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell + this->p_hamilt = new hamilt::HamiltLIP(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell #ifdef __EXX , *this->exx_lip #endif diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h index 6684281689..56276f3eb6 100644 --- a/source/module_esolver/esolver_ks_lcaopw.h +++ b/source/module_esolver/esolver_ks_lcaopw.h @@ -33,7 +33,7 @@ namespace ModuleESolver const int iter, const double ethr) override; - virtual void allocate_hamilt() override; + virtual void allocate_hamilt(const UnitCell& ucell) override; virtual void deallocate_hamilt() override; #ifdef __EXX diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index b70c53eec0..67e4d33dd1 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -115,9 +115,9 @@ ESolver_KS_PW::~ESolver_KS_PW() } template -void ESolver_KS_PW::allocate_hamilt() +void ESolver_KS_PW::allocate_hamilt(const UnitCell& ucell) { - this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell); + this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell); } template void ESolver_KS_PW::deallocate_hamilt() @@ -202,10 +202,10 @@ void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_p } //! init pseudopotential - this->ppcell.init(ucell.ntype, &this->sf, this->pw_wfc); + this->ppcell.init(ucell,&this->sf, this->pw_wfc); //! initalize local pseudopotential - this->ppcell.init_vloc(this->pw_rhod); + this->ppcell.init_vloc(ucell,this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); //! Initalize non-local pseudopotential @@ -279,7 +279,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->deallocate_hamilt(); // allocate HamiltPW - this->allocate_hamilt(); + this->allocate_hamilt(ucell); //---------------------------------------------------------- // about vdw, jiyy add vdwd3 and linpz add vdwd2 @@ -621,7 +621,8 @@ void ESolver_KS_PW::cal_force(UnitCell& ucell, ModuleBase::matrix& fo : reinterpret_cast, Device>*>(this->kspw_psi); // Calculate forces - ff.cal_force(force, + ff.cal_force(ucell, + force, *this->pelec, this->pw_rhod, &ucell.symm, diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index da5a7dcc22..3ad3aeb09c 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -46,7 +46,7 @@ class ESolver_KS_PW : public ESolver_KS virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override; - virtual void allocate_hamilt(); + virtual void allocate_hamilt(const UnitCell& ucell); virtual void deallocate_hamilt(); //! hide the psi in ESolver_KS for tmp use diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 74dd4fa98d..32d4a90b2e 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -112,7 +112,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); // initialize local pseudopotential - this->locpp.init_vloc(pw_rho); + this->locpp.init_vloc(ucell,pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); @@ -538,7 +538,7 @@ double ESolver_OF::cal_energy() void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force) { Forces ff(ucell.nat); - ff.cal_force(force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp); + ff.cal_force(ucell,force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp); } /** diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 76b510a36a..fc18e8e302 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -100,6 +100,7 @@ void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) this->pw_wfc, &this->kv, &this->ppcell, + &ucell, PARAM.globalv.npol, &this->stoche.emin_sto, &this->stoche.emax_sto); @@ -166,7 +167,8 @@ void ESolver_SDFT_PW::hamilt2density_single(UnitCell& ucell, int iste hsolver::DiagoIterAssist::PW_DIAG_THR, hsolver::DiagoIterAssist::need_subspace); - hsolver_pw_sdft_obj.solve(this->p_hamilt, + hsolver_pw_sdft_obj.solve(ucell, + this->p_hamilt, this->kspw_psi[0], this->psi[0], this->pelec, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index b214cfe8a0..389119a17a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -869,11 +869,11 @@ void Force_Stress_LCAO::calForcePwPart(const UnitCell& ucell, // local pseudopotential force: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - f_pw.cal_force_loc(fvl_dvl, rhopw, nlpp.vloc, chr); + f_pw.cal_force_loc(ucell,fvl_dvl, rhopw, nlpp.vloc, chr); //-------------------------------------------------------- // ewald force: use plane wave only. //-------------------------------------------------------- - f_pw.cal_force_ew(fewalds, rhopw, &sf); // remain problem + f_pw.cal_force_ew(ucell,fewalds, rhopw, &sf); // remain problem //-------------------------------------------------------- // force due to core correlation. @@ -1009,17 +1009,17 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, // local pseudopotential stress: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - sc_pw.stress_loc(sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr); + sc_pw.stress_loc(ucell,sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr); //-------------------------------------------------------- // hartree term //-------------------------------------------------------- - sc_pw.stress_har(sigmahar, rhopw, 0, chr); + sc_pw.stress_har(ucell,sigmahar, rhopw, 0, chr); //-------------------------------------------------------- // ewald stress: use plane wave only. //-------------------------------------------------------- - sc_pw.stress_ewa(sigmaewa, rhopw, 0); // remain problem + sc_pw.stress_ewa(ucell,sigmaewa, rhopw, 0); // remain problem //-------------------------------------------------------- // stress due to core correlation. @@ -1034,7 +1034,7 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, sigmaxc(i, i) = -etxc / ucell.omega; } // Exchange-correlation for PBE - sc_pw.stress_gga(sigmaxc, rhopw, chr); + sc_pw.stress_gga(ucell,sigmaxc, rhopw, chr); return; } 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 1f28939c09..d0f791d34c 100644 --- a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp @@ -63,20 +63,20 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, } // hartree contribution - stress_har(sigmahar, this->rhopw, true, pelec->charge); + stress_har(ucell,sigmahar, this->rhopw, true, pelec->charge); // ewald contribution - stress_ewa(sigmaewa, this->rhopw, true); + stress_ewa(ucell,sigmaewa, this->rhopw, true); // xc contribution: add gradient corrections(non diagonal) for (int i = 0; i < 3; i++) { sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega; } - stress_gga(sigmaxc, this->rhopw, pelec->charge); + stress_gga(ucell,sigmaxc, this->rhopw, pelec->charge); // local contribution - stress_loc(sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge); + 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); diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp index 98f4198036..486b866bf1 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp @@ -17,7 +17,8 @@ pseudopot_cell_vl::~pseudopot_cell_vl() delete[] zp; } -void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) +void pseudopot_cell_vl::init_vloc(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis) { if(PARAM.inp.use_paw) { return; } @@ -30,11 +31,11 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) double *vloc1d = new double[rho_basis->ngg]; ModuleBase::GlobalFunc::ZEROS(vloc1d, rho_basis->ngg); - this->allocate(rho_basis->ngg); + this->allocate(ucell,rho_basis->ngg); - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - const Atom* atom = &GlobalC::ucell.atoms[it]; + const Atom* atom = &ucell.atoms[it]; ModuleBase::GlobalFunc::ZEROS(vloc1d, rho_basis->ngg); @@ -42,7 +43,7 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) // compute V_loc(G) for a given type of atom if(atom->coulomb_potential) { - this->vloc_coulomb(this->zp[it], vloc1d, rho_basis); + this->vloc_coulomb(ucell,this->zp[it], vloc1d, rho_basis); } else if(numeric[it]==true) { @@ -53,6 +54,7 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) atom->ncpp.vloc_at.data(), // local potential in real space radial form. this->zp[it], vloc1d, + ucell, rho_basis); } else @@ -69,26 +71,27 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis) delete[] vloc1d; - this->print_vloc(rho_basis); + this->print_vloc(ucell,rho_basis); ModuleBase::timer::tick("ppcell_vl","init_vloc"); return; } -void pseudopot_cell_vl::allocate(const int ngg) +void pseudopot_cell_vl::allocate(const UnitCell& ucell, + const int ngg) { if(PARAM.inp.test_pp>0) { ModuleBase::TITLE("pseudopot_cell_vl","allocate"); } if(PARAM.inp.use_paw) { return; } - this->vloc.create(GlobalC::ucell.ntype, ngg); + this->vloc.create(ucell.ntype, ngg); delete[] numeric; - this->numeric = new bool[GlobalC::ucell.ntype]; - ModuleBase::GlobalFunc::ZEROS(numeric, GlobalC::ucell.ntype); + this->numeric = new bool[ucell.ntype]; + ModuleBase::GlobalFunc::ZEROS(numeric, ucell.ntype); - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { this->numeric[it] = true; } @@ -104,7 +107,10 @@ void pseudopot_cell_vl::allocate(const int ngg) return; } -void pseudopot_cell_vl::vloc_coulomb(const double& zp_in, double* vloc_1d, const ModulePW::PW_Basis* rho_basis) const +void pseudopot_cell_vl::vloc_coulomb(const UnitCell& ucell, + const double& zp_in, + double* vloc_1d, + const ModulePW::PW_Basis* rho_basis) const { int igl0 = 0; // start from |G|=0 or not. @@ -117,14 +123,14 @@ void pseudopot_cell_vl::vloc_coulomb(const double& zp_in, double* vloc_1d, const { igl0 = 0; } - const double d_fpi_omega = ModuleBase::FOUR_PI / GlobalC::ucell.omega; // mohan add 2008-06-04 + const double d_fpi_omega = ModuleBase::FOUR_PI / ucell.omega; // mohan add 2008-06-04 double fac = -zp_in * ModuleBase::e2 * d_fpi_omega; #ifdef _OPENMP #pragma omp for #endif for (int ig = igl0; ig < rho_basis->ngg; ig++) { - double gx2 = rho_basis->gg_uniq[ig] * GlobalC::ucell.tpiba2; + double gx2 = rho_basis->gg_uniq[ig] * ucell.tpiba2; vloc_1d[ig] = fac / gx2; } return; @@ -145,6 +151,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh, const double* vloc_at, const double& zp_in, double* vloc_1d, + const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis) const { //---------------------------------------------------------------- @@ -173,7 +180,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh, /* for(ir=0; irngg;ig++) { - double gx2= rho_basis->gg_uniq[ig] * GlobalC::ucell.tpiba2; + double gx2= rho_basis->gg_uniq[ig] * ucell.tpiba2; double gx = std::sqrt(gx2); for (int ir = 0;ir < msh;ir++) { @@ -235,7 +242,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh, vloc_1d[ig] -= fac * ModuleBase::libm::exp(- gx2 * 0.25)/ gx2; } // enddo - const double d_fpi_omega = ModuleBase::FOUR_PI/GlobalC::ucell.omega;//mohan add 2008-06-04 + const double d_fpi_omega = ModuleBase::FOUR_PI/ucell.omega;//mohan add 2008-06-04 #ifdef _OPENMP #pragma omp for #endif @@ -253,21 +260,22 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh, return; } // end subroutine vloc_of_g -void pseudopot_cell_vl::print_vloc(const ModulePW::PW_Basis* rho_basis) const +void pseudopot_cell_vl::print_vloc(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis) const { if(GlobalV::MY_RANK!=0) { return; //mohan fix bug 2011-10-13 } bool check_vl = PARAM.inp.out_element_info; if(check_vl) { - for(int it=0; itngg;ig++) { - ofs_vg << std::setw(15) << rho_basis->gg_uniq [ig] * GlobalC::ucell.tpiba2 + ofs_vg << std::setw(15) << rho_basis->gg_uniq [ig] * ucell.tpiba2 << std::setw(15) << this->vloc(it, ig) << std::endl; } ofs_vg.close(); diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h index a5f91b0f81..5d439216bb 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h @@ -5,6 +5,7 @@ #include "module_base/global_variable.h" #include "module_base/matrix.h" #include "module_basis/module_pw/pw_basis.h" +#include "module_cell/unitcell.h" class pseudopot_cell_vl { @@ -19,7 +20,8 @@ class pseudopot_cell_vl * @param rho_basis pw basis * @return this->vloc */ - void init_vloc(const ModulePW::PW_Basis* rho_basis); + void init_vloc(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis); ModuleBase::matrix vloc; //(ntype,ngl),the local potential for each atom type(ntype,ngl) bool *numeric; //[ntype], =true @@ -28,12 +30,16 @@ class pseudopot_cell_vl double *zp; // (npsx),the charge of the pseudopotential - void allocate(const int ngg); + void allocate(const UnitCell& ucell, + const int ngg); /** * @brief compute the coulomb potential in reciprocal space * v(g) = -\frac{4pi}{V} * zp*e^2 / G^2 */ - void vloc_coulomb(const double& zp, double* vloc_1d, const ModulePW::PW_Basis* rho_basis) const; + void vloc_coulomb(const UnitCell& ucell, + const double& zp, + double* vloc_1d, + const ModulePW::PW_Basis* rho_basis) const; // generate vloc for a particular atom type. void vloc_of_g(const int& msh, const double* rab, @@ -41,9 +47,11 @@ class pseudopot_cell_vl const double* vloc_at, const double& zp, double* vloc, + const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis) const; - void print_vloc(const ModulePW::PW_Basis* rho_basis) const; + void print_vloc(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis) const; }; #endif // VL_IN_PW diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp index 161f72e0f5..856ca038f7 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp @@ -58,7 +58,8 @@ void pseudopot_cell_vnl::initgradq_vnl(const UnitCell &cell) } -void pseudopot_cell_vnl::getgradq_vnl(const int ik) +void pseudopot_cell_vnl::getgradq_vnl(const UnitCell& ucell, + const int ik) { if(PARAM.inp.test_pp) ModuleBase::TITLE("pseudopot_cell_vnl","getvnl"); ModuleBase::timer::tick("pp_cell_vnl","getvnl"); @@ -91,11 +92,11 @@ void pseudopot_cell_vnl::getgradq_vnl(const int ik) ModuleBase::YlmReal::grad_Ylm_Real(x1, npw, gk, ylm, dylm[0], dylm[1], dylm[2]); int jkb = 0; - for(int it = 0;it < GlobalC::ucell.ntype;it++) + for(int it = 0;it < ucell.ntype;it++) { // calculate beta in G-space using an interpolation table - const int nbeta = GlobalC::ucell.atoms[it].ncpp.nbeta; - const int nh = GlobalC::ucell.atoms[it].ncpp.nh; + const int nbeta = ucell.atoms[it].ncpp.nbeta; + const int nh = ucell.atoms[it].ncpp.nh; int nb0 = -1; for( int ih = 0; ih < nh; ++ih) { @@ -104,7 +105,7 @@ void pseudopot_cell_vnl::getgradq_vnl(const int ik) { for (int ig = 0;ig < npw;++ig) { - const double gnorm = gk[ig].norm() * GlobalC::ucell.tpiba; + const double gnorm = gk[ig].norm() * ucell.tpiba; vq [ig] = ModuleBase::PolyInt::Polynomial_Interpolation( this->tab, it, nb, PARAM.globalv.nqx, PARAM.globalv.dq, gnorm ); dvq[ig] =ModuleBase::PolyInt::Polynomial_Interpolation( @@ -146,7 +147,7 @@ void pseudopot_cell_vnl::getgradq_vnl(const int ik) // vkb1 contains all betas including angular part for type nt // now add the structure factor and factor (-i)^l - for (int ia=0; ia *sk = this->psf->get_sk(ik, it, ia, this->wfcpw); @@ -161,7 +162,7 @@ void pseudopot_cell_vnl::getgradq_vnl(const int ik) { std::complex skig = sk[ig]; pvkb[ig] = tmpvkb(ih, ig) * skig * pref; - // std::complex dskig = ModuleBase::NEG_IMAG_UNIT * (GlobalC::ucell.atoms[it].tau[ia][id] * this->wfcpw->lat0) * skig; + // std::complex dskig = ModuleBase::NEG_IMAG_UNIT * (ucell.atoms[it].tau[ia][id] * this->wfcpw->lat0) * skig; // pgvkb[ig] = tmpgradvkb(id, ih, ig) * skig * pref + tmpvkb(ih, ig) * dskig * pref; // The second term will be eliminate when doing Dij or we can say (\nabla_q+\nabla_q')S(q'-q) = 0 pgvkb[ig] = tmpgradvkb(id, ih, ig) * skig * pref; diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp index dee8205a42..b1307b8523 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp @@ -82,15 +82,16 @@ void pseudopot_cell_vnl::release_memory() // setup lmaxkb, nhm, nkb, lmaxq // allocate vkb, PARAM.globalv.nqx, tab, tab_at //----------------------------------- -void pseudopot_cell_vnl::init(const int ntype, +void pseudopot_cell_vnl::init(const UnitCell& ucell, Structure_Factor* psf_in, const ModulePW::PW_Basis_K* wfc_basis, const bool allocate_vkb) { - if (PARAM.inp.use_paw) { + if (PARAM.inp.use_paw) + { return; -} - + } + const int ntype = ucell.ntype; ModuleBase::TITLE("pseudopot_cell_vnl", "init"); ModuleBase::timer::tick("ppcell_vnl", "init"); @@ -106,12 +107,12 @@ void pseudopot_cell_vnl::init(const int ntype, this->lmaxkb = -1; for (it = 0; it < ntype; it++) { - GlobalV::ofs_running << " " << GlobalC::ucell.atoms[it].label << " non-local projectors:" << std::endl; - for (int ibeta = 0; ibeta < GlobalC::ucell.atoms[it].ncpp.nbeta; ibeta++) + GlobalV::ofs_running << " " << ucell.atoms[it].label << " non-local projectors:" << std::endl; + for (int ibeta = 0; ibeta < ucell.atoms[it].ncpp.nbeta; ibeta++) { - GlobalV::ofs_running << " projector " << ibeta + 1 << " L=" << GlobalC::ucell.atoms[it].ncpp.lll[ibeta] + GlobalV::ofs_running << " projector " << ibeta + 1 << " L=" << ucell.atoms[it].ncpp.lll[ibeta] << std::endl; - this->lmaxkb = std::max(this->lmaxkb, GlobalC::ucell.atoms[it].ncpp.lll[ibeta]); + this->lmaxkb = std::max(this->lmaxkb, ucell.atoms[it].ncpp.lll[ibeta]); } } @@ -126,9 +127,9 @@ void pseudopot_cell_vnl::init(const int ntype, int nwfcm = 0; for (it = 0; it < ntype; it++) { - this->nhm = std::max(nhm, GlobalC::ucell.atoms[it].ncpp.nh); - this->nbetam = std::max(nbetam, GlobalC::ucell.atoms[it].ncpp.nbeta); - nwfcm = std::max(nwfcm, GlobalC::ucell.atoms[it].ncpp.nchi); + this->nhm = std::max(nhm, ucell.atoms[it].ncpp.nh); + this->nbetam = std::max(nbetam, ucell.atoms[it].ncpp.nbeta); + nwfcm = std::max(nwfcm, ucell.atoms[it].ncpp.nchi); } //---------------------------------------------------------- @@ -138,7 +139,7 @@ void pseudopot_cell_vnl::init(const int ntype, this->nkb = 0; for (it = 0; it < ntype; it++) { - this->nkb += GlobalC::ucell.atoms[it].ncpp.nh * GlobalC::ucell.atoms[it].na; + this->nkb += ucell.atoms[it].ncpp.nh * ucell.atoms[it].na; } ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "TOTAL NUMBER OF NONLOCAL PROJECTORS", nkb); @@ -149,28 +150,28 @@ void pseudopot_cell_vnl::init(const int ntype, this->nhtol.create(ntype, this->nhm); this->nhtolm.create(ntype, this->nhm); this->nhtoj.create(ntype, this->nhm); - this->deeq.create(PARAM.inp.nspin, GlobalC::ucell.nat, this->nhm, this->nhm); - this->deeq_nc.create(PARAM.inp.nspin, GlobalC::ucell.nat, this->nhm, this->nhm); + this->deeq.create(PARAM.inp.nspin, ucell.nat, this->nhm, this->nhm); + this->deeq_nc.create(PARAM.inp.nspin, ucell.nat, this->nhm, this->nhm); this->qq_nt.create(ntype, this->nhm, this->nhm); this->qq_so.create(ntype, 4, this->nhm, this->nhm); if (PARAM.inp.device == "gpu") { if (PARAM.inp.precision == "single") { - resmem_sd_op()(gpu_ctx, s_deeq, PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm); + resmem_sd_op()(gpu_ctx, s_deeq, PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm); resmem_sd_op()(gpu_ctx, s_nhtol, ntype * this->nhm); resmem_sd_op()(gpu_ctx, s_nhtolm, ntype * this->nhm); resmem_sd_op()(gpu_ctx, s_indv, ntype * this->nhm); resmem_sd_op()(gpu_ctx, s_qq_nt, ntype * this->nhm * this->nhm); - resmem_cd_op()(gpu_ctx, c_deeq_nc, PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm); + resmem_cd_op()(gpu_ctx, c_deeq_nc, PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm); resmem_cd_op()(gpu_ctx, c_qq_so, ntype * 4 * this->nhm * this->nhm); } else { - resmem_zd_op()(gpu_ctx, z_deeq_nc, PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm); + resmem_zd_op()(gpu_ctx, z_deeq_nc, PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm); resmem_zd_op()(gpu_ctx, z_qq_so, ntype * 4 * this->nhm * this->nhm); } - resmem_dd_op()(gpu_ctx, d_deeq, PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm); + resmem_dd_op()(gpu_ctx, d_deeq, PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm); resmem_dd_op()(gpu_ctx, d_indv, ntype * this->nhm); resmem_dd_op()(gpu_ctx, d_nhtol, ntype * this->nhm); resmem_dd_op()(gpu_ctx, d_nhtolm, ntype * this->nhm); @@ -182,7 +183,7 @@ void pseudopot_cell_vnl::init(const int ntype, { resmem_sh_op()(cpu_ctx, s_deeq, - PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm, + PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm, "VNL::s_deeq"); resmem_sh_op()(cpu_ctx, s_nhtol, ntype * this->nhm, "VNL::s_nhtol"); resmem_sh_op()(cpu_ctx, s_nhtolm, ntype * this->nhm, "VNL::s_nhtolm"); @@ -190,7 +191,7 @@ void pseudopot_cell_vnl::init(const int ntype, resmem_sh_op()(cpu_ctx, s_qq_nt, ntype * this->nhm * this->nhm, "VNL::s_qq_nt"); resmem_ch_op()(cpu_ctx, c_deeq_nc, - PARAM.inp.nspin * GlobalC::ucell.nat * this->nhm * this->nhm, + PARAM.inp.nspin * ucell.nat * this->nhm * this->nhm, "VNL::c_deeq_nc"); resmem_ch_op()(cpu_ctx, c_qq_so, ntype * 4 * this->nhm * this->nhm, "VNL::c_qq_so"); } @@ -210,7 +211,7 @@ void pseudopot_cell_vnl::init(const int ntype, this->dvan_so.create(PARAM.inp.nspin, ntype, this->nhm, this->nhm); this->ijtoh.create(ntype, this->nhm, this->nhm); - this->qq_at.create(GlobalC::ucell.nat, this->nhm, this->nhm); + this->qq_at.create(ucell.nat, this->nhm, this->nhm); } else { @@ -300,7 +301,7 @@ void pseudopot_cell_vnl::init(const int ntype, // Calculates beta functions (Kleinman-Bylander projectors), // with structure factor, for all atoms, in reciprocal space //---------------------------------------------------------- -void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in) const +void pseudopot_cell_vnl::getvnl(const int& ik, const UnitCell& ucell, ModuleBase::ComplexMatrix& vkb_in) const { if (PARAM.inp.use_paw) { return; @@ -338,15 +339,15 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in using resmem_complex_op = base_device::memory::resize_memory_op, Device>; using delmem_complex_op = base_device::memory::delete_memory_op, Device>; std::complex* sk = nullptr; - resmem_complex_op()(ctx, sk, GlobalC::ucell.nat * npw, "VNL::sk"); + resmem_complex_op()(ctx, sk, ucell.nat * npw, "VNL::sk"); this->psf->get_sk(ctx, ik, this->wfcpw, sk); int jkb = 0, iat = 0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { // calculate beta in G-space using an interpolation table - const int nbeta = GlobalC::ucell.atoms[it].ncpp.nbeta; - const int nh = GlobalC::ucell.atoms[it].ncpp.nh; + const int nbeta = ucell.atoms[it].ncpp.nbeta; + const int nh = ucell.atoms[it].ncpp.nh; if (PARAM.inp.test_pp > 1) { ModuleBase::GlobalFunc::OUT("nbeta", nbeta); @@ -359,7 +360,7 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in } for (int ig = 0; ig < npw; ig++) { - const double gnorm = gk[ig].norm() * GlobalC::ucell.tpiba; + const double gnorm = gk[ig].norm() * ucell.tpiba; vq[ig] = ModuleBase::PolyInt::Polynomial_Interpolation(this->tab, it, @@ -385,7 +386,7 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in // vkb1 contains all betas including angular part for type nt // now add the structure factor and factor (-i)^l - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int ia = 0; ia < ucell.atoms[it].na; ia++) { for (int ih = 0; ih < nh; ih++) { @@ -410,7 +411,10 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in } // end subroutine getvnl template -void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex* vkb_in) const +void pseudopot_cell_vnl::getvnl(Device* ctx, + const UnitCell& ucell, + const int& ik, + std::complex* vkb_in) const { if (PARAM.inp.use_paw) { return; @@ -440,13 +444,13 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex const int x1 = (lmaxkb + 1) * (lmaxkb + 1); const int npw = this->wfcpw->npwk[ik]; - int *atom_nh = nullptr, *atom_na = nullptr, *atom_nb = nullptr, *h_atom_nh = new int[GlobalC::ucell.ntype], - *h_atom_na = new int[GlobalC::ucell.ntype], *h_atom_nb = new int[GlobalC::ucell.ntype]; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + int *atom_nh = nullptr, *atom_na = nullptr, *atom_nb = nullptr, *h_atom_nh = new int[ucell.ntype], + *h_atom_na = new int[ucell.ntype], *h_atom_nb = new int[ucell.ntype]; + for (int it = 0; it < ucell.ntype; it++) { - h_atom_nb[it] = GlobalC::ucell.atoms[it].ncpp.nbeta; - h_atom_nh[it] = GlobalC::ucell.atoms[it].ncpp.nh; - h_atom_na[it] = GlobalC::ucell.atoms[it].na; + h_atom_nb[it] = ucell.atoms[it].ncpp.nbeta; + h_atom_nh[it] = ucell.atoms[it].ncpp.nh; + h_atom_na[it] = ucell.atoms[it].na; } // When the internal memory is large enough, it is better to make vkb1 be the number of pseudopot_cell_vnl. // We only need to initialize it once as long as the cell is unchanged. @@ -466,12 +470,12 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex } if (PARAM.inp.device == "gpu") { - resmem_int_op()(ctx, atom_nh, GlobalC::ucell.ntype); - resmem_int_op()(ctx, atom_nb, GlobalC::ucell.ntype); - resmem_int_op()(ctx, atom_na, GlobalC::ucell.ntype); - syncmem_int_op()(ctx, cpu_ctx, atom_nh, h_atom_nh, GlobalC::ucell.ntype); - syncmem_int_op()(ctx, cpu_ctx, atom_nb, h_atom_nb, GlobalC::ucell.ntype); - syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, GlobalC::ucell.ntype); + resmem_int_op()(ctx, atom_nh, ucell.ntype); + resmem_int_op()(ctx, atom_nb, ucell.ntype); + resmem_int_op()(ctx, atom_na, ucell.ntype); + syncmem_int_op()(ctx, cpu_ctx, atom_nh, h_atom_nh, ucell.ntype); + syncmem_int_op()(ctx, cpu_ctx, atom_nb, h_atom_nb, ucell.ntype); + syncmem_int_op()(ctx, cpu_ctx, atom_na, h_atom_na, ucell.ntype); resmem_var_op()(ctx, gk, npw * 3); castmem_var_h2d_op()(ctx, cpu_ctx, gk, reinterpret_cast(_gk), npw * 3); @@ -495,11 +499,11 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex ModuleBase::YlmReal::Ylm_Real(ctx, x1, npw, gk, ylm); std::complex* sk = nullptr; - resmem_complex_op()(ctx, sk, GlobalC::ucell.nat * npw); + resmem_complex_op()(ctx, sk, ucell.nat * npw); this->psf->get_sk(ctx, ik, this->wfcpw, sk); cal_vnl_op()(ctx, - GlobalC::ucell.ntype, + ucell.ntype, npw, this->wfcpw->npwk_max, this->nhm, @@ -509,7 +513,7 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex atom_nb, atom_nh, static_cast(PARAM.globalv.dq), - static_cast(GlobalC::ucell.tpiba), + static_cast(ucell.tpiba), static_cast>(ModuleBase::NEG_IMAG_UNIT), gk, ylm, @@ -584,7 +588,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_ this->dvan.zero_out(); this->dvan_so.zero_out(); // added by zhengdy-soc delete[] indv_ijkb0; - this->indv_ijkb0 = new int[GlobalC::ucell.nat]; + this->indv_ijkb0 = new int[cell.nat]; int ijkb0 = 0; for (int it = 0; it < cell.ntype; it++) { @@ -1219,12 +1223,12 @@ double pseudopot_cell_vnl::CG(int l1, int m1, int l2, int m2, int L, int M) // p // MGT.init_Gaunt( lmaxkb + 2 ); // int jkb = 0; -// for(int it = 0;it < GlobalC::ucell.ntype;it++) +// for(int it = 0;it < ucell.ntype;it++) // { // if(PARAM.inp.test_pp>1) ModuleBase::GlobalFunc::OUT("it",it); // // calculate beta in G-space using an interpolation table -// const int nbeta = GlobalC::ucell.atoms[it].ncpp.nbeta; -// const int nh = GlobalC::ucell.atoms[it].ncpp.nh; +// const int nbeta = ucell.atoms[it].ncpp.nbeta; +// const int nh = ucell.atoms[it].ncpp.nh; // if(PARAM.inp.test_pp>1) ModuleBase::GlobalFunc::OUT("nbeta",nbeta); @@ -1244,7 +1248,7 @@ double pseudopot_cell_vnl::CG(int l1, int m1, int l2, int m2, int L, int M) // p // { // for (ig = 0;ig < npw;ig++) // { -// const double gnorm = gk[ig].norm() * GlobalC::ucell.tpiba; +// const double gnorm = gk[ig].norm() * ucell.tpiba; // vq [ig] = ModuleBase::PolyInt::Polynomial_Interpolation( // this->tab_alpha, it, nb, L, PARAM.globalv.nqx, PARAM.globalv.dq, gnorm); @@ -1266,7 +1270,7 @@ double pseudopot_cell_vnl::CG(int l1, int m1, int l2, int m2, int L, int M) // p // } // } // end nbeta -// for (ia=0; ia *sk = this->psf->get_sk(ik, it, ia,this->wfcpw); // for (ih = 0;ih < nh;ih++) @@ -1291,20 +1295,20 @@ double pseudopot_cell_vnl::CG(int l1, int m1, int l2, int m2, int L, int M) // p // } #endif -void pseudopot_cell_vnl::init_vnl_alpha() // pengfei Li 2018-3-23 +void pseudopot_cell_vnl::init_vnl_alpha(const UnitCell& ucell) // pengfei Li 2018-3-23 { if (PARAM.inp.test_pp) { ModuleBase::TITLE("pseudopot_cell_vnl", "init_vnl_alpha"); } ModuleBase::timer::tick("ppcell_vnl", "init_vnl_alpha"); - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { int BetaIndex = 0; - // const int Nprojectors = GlobalC::ucell.atoms[it].nh; - for (int ib = 0; ib < GlobalC::ucell.atoms[it].ncpp.nbeta; ib++) + // const int Nprojectors = ucell.atoms[it].nh; + for (int ib = 0; ib < ucell.atoms[it].ncpp.nbeta; ib++) { - const int l = GlobalC::ucell.atoms[it].ncpp.lll[ib]; + const int l = ucell.atoms[it].ncpp.lll[ib]; for (int m = 0; m < 2 * l + 1; m++) { this->nhtol(it, BetaIndex) = l; @@ -1318,14 +1322,14 @@ void pseudopot_cell_vnl::init_vnl_alpha() // pengfei Li 2018-3-23 // max number of beta functions const int nbrx = 10; - const double pref = ModuleBase::FOUR_PI / sqrt(GlobalC::ucell.omega); - this->tab_alpha.create(GlobalC::ucell.ntype, nbrx, lmaxkb + 2, PARAM.globalv.nqx); + const double pref = ModuleBase::FOUR_PI / sqrt(ucell.omega); + this->tab_alpha.create(ucell.ntype, nbrx, lmaxkb + 2, PARAM.globalv.nqx); this->tab_alpha.zero_out(); GlobalV::ofs_running << "\n Init Non-Local PseudoPotential table( including L index) : "; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - const int nbeta = GlobalC::ucell.atoms[it].ncpp.nbeta; - int kkbeta = GlobalC::ucell.atoms[it].ncpp.kkbeta; + const int nbeta = ucell.atoms[it].ncpp.nbeta; + int kkbeta = ucell.atoms[it].ncpp.kkbeta; // mohan modify 2008-3-31 // mohan add kkbeta>0 2009-2-27 @@ -1344,15 +1348,15 @@ void pseudopot_cell_vnl::init_vnl_alpha() // pengfei Li 2018-3-23 for (int iq = 0; iq < PARAM.globalv.nqx; iq++) { const double q = iq * PARAM.globalv.dq; - ModuleBase::Sphbes::Spherical_Bessel(kkbeta, GlobalC::ucell.atoms[it].ncpp.r.data(), q, L, jl); + ModuleBase::Sphbes::Spherical_Bessel(kkbeta, ucell.atoms[it].ncpp.r.data(), q, L, jl); for (int ir = 0; ir < kkbeta; ir++) { - aux[ir] = GlobalC::ucell.atoms[it].ncpp.betar(ib, ir) * jl[ir] - * GlobalC::ucell.atoms[it].ncpp.r[ir] * GlobalC::ucell.atoms[it].ncpp.r[ir]; + aux[ir] = ucell.atoms[it].ncpp.betar(ib, ir) * jl[ir] + * ucell.atoms[it].ncpp.r[ir] * ucell.atoms[it].ncpp.r[ir]; } double vqint; - ModuleBase::Integral::Simpson_Integral(kkbeta, aux, GlobalC::ucell.atoms[it].ncpp.rab.data(), vqint); + ModuleBase::Integral::Simpson_Integral(kkbeta, aux, ucell.atoms[it].ncpp.rab.data(), vqint); this->tab_alpha(it, ib, L, iq) = vqint * pref; } } @@ -1827,16 +1831,20 @@ std::complex* pseudopot_cell_vnl::get_qq_so_data() const } template void pseudopot_cell_vnl::getvnl(base_device::DEVICE_CPU*, + const UnitCell&, int const&, std::complex*) const; template void pseudopot_cell_vnl::getvnl(base_device::DEVICE_CPU*, + const UnitCell&, int const&, std::complex*) const; #if defined(__CUDA) || defined(__ROCM) template void pseudopot_cell_vnl::getvnl(base_device::DEVICE_GPU*, + const UnitCell&, int const&, std::complex*) const; template void pseudopot_cell_vnl::getvnl(base_device::DEVICE_GPU*, + const UnitCell&, int const&, std::complex*) const; #endif diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h index d3fa7d2164..e647675e87 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h @@ -24,10 +24,10 @@ class pseudopot_cell_vnl : public pseudopot_cell_vl public: pseudopot_cell_vnl(); ~pseudopot_cell_vnl(); - void init(const int ntype, + void init(const UnitCell& cell, Structure_Factor* psf_in, const ModulePW::PW_Basis_K* wfc_basis = nullptr, - const bool allocate_vkb = 1); + const bool allocate_vkb = true); double cell_factor; // LiuXh add 20180619 @@ -38,17 +38,23 @@ class pseudopot_cell_vnl : public pseudopot_cell_vl void init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_basis); template - void getvnl(Device* ctx, const int& ik, std::complex* vkb_in) const; + void getvnl(Device* ctx, + const UnitCell& ucell, + const int& ik, + std::complex* vkb_in) const; - void getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in) const; + void getvnl(const int& ik, + const UnitCell& ucell, + ModuleBase::ComplexMatrix& vkb_in) const; // void getvnl_alpha(const int &ik); - void init_vnl_alpha(void); + void init_vnl_alpha(const UnitCell& cell); void initgradq_vnl(const UnitCell& cell); - void getgradq_vnl(const int ik); + void getgradq_vnl(const UnitCell& ucell, + const int ik); //=============================================================== // MEMBER VARIABLES : diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index 0ffce5313a..f9c6a63556 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -24,7 +24,8 @@ #endif template -void Forces::cal_force(ModuleBase::matrix& force, +void Forces::cal_force(const UnitCell& ucell, + ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, @@ -54,7 +55,7 @@ void Forces::cal_force(ModuleBase::matrix& force, // For PAW, calculated together in paw_cell.calculate_force if (!PARAM.inp.use_paw) { - this->cal_force_loc(forcelc, rho_basis, locpp->vloc, chr); + this->cal_force_loc(ucell,forcelc, rho_basis, locpp->vloc, chr); } else { @@ -62,7 +63,7 @@ void Forces::cal_force(ModuleBase::matrix& force, } // Ewald - this->cal_force_ew(forceion, rho_basis, p_sf); + this->cal_force_ew(ucell,forceion, rho_basis, p_sf); // Force due to nonlocal part of pseudopotential if (wfc_basis != nullptr) @@ -70,11 +71,11 @@ void Forces::cal_force(ModuleBase::matrix& force, if (!PARAM.inp.use_paw) { this->npwx = wfc_basis->npwk_max; - Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, *p_nlpp, GlobalC::ucell, psi_in); + Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, *p_nlpp, ucell, psi_in); if (PARAM.globalv.use_uspp) { - this->cal_force_us(forcenl, rho_basis, *p_nlpp, elec, GlobalC::ucell); + this->cal_force_us(forcenl, rho_basis, *p_nlpp, elec, ucell); } } else @@ -119,7 +120,7 @@ void Forces::cal_force(ModuleBase::matrix& force, wfc_basis->get_ig2iy(ik).data(), wfc_basis->get_ig2iz(ik).data(), (const double**)kpg, - GlobalC::ucell.tpiba, + ucell.tpiba, (const double**)gcar); delete[] kpt; @@ -172,7 +173,7 @@ void Forces::cal_force(ModuleBase::matrix& force, // For PAW, calculated together in paw_cell.calculate_force if (!PARAM.inp.use_paw) { - this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, GlobalC::ucell); + this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell); } else { @@ -182,7 +183,7 @@ void Forces::cal_force(ModuleBase::matrix& force, ModuleBase::matrix stress_vdw_pw; //.create(3,3); ModuleBase::matrix force_vdw; force_vdw.create(nat, 3); - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); + auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); if (vdw_solver != nullptr) { const std::vector>& force_vdw_temp = vdw_solver->get_force(); @@ -194,7 +195,7 @@ void Forces::cal_force(ModuleBase::matrix& force, } if (PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "VDW FORCE (Ry/Bohr)", force_vdw); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "VDW FORCE (Ry/Bohr)", force_vdw); } } @@ -202,10 +203,10 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.efield_flag) { force_e.create(this->nat, 3); - elecstate::Efield::compute_force(GlobalC::ucell, force_e); + elecstate::Efield::compute_force(ucell, force_e); if (PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD FORCE (Ry/Bohr)", force_e); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD FORCE (Ry/Bohr)", force_e); } } @@ -213,10 +214,10 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.gate_flag) { force_gate.create(this->nat, 3); - elecstate::Gatefield::compute_force(GlobalC::ucell, force_gate); + elecstate::Gatefield::compute_force(ucell, force_gate); if (PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "GATEFIELD FORCE (Ry/Bohr)", force_gate); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "GATEFIELD FORCE (Ry/Bohr)", force_gate); } } @@ -224,10 +225,10 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.imp_sol) { forcesol.create(this->nat, 3); - GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rho_basis, locpp->vloc, forcesol); + GlobalC::solvent_model.cal_force_sol(ucell, rho_basis, locpp->vloc, forcesol); if (PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); } } @@ -253,13 +254,13 @@ void Forces::cal_force(ModuleBase::matrix& force, ModuleBase::matrix v_xc, v_effective; v_effective.create(PARAM.inp.nspin, rho_basis->nrxx); v_effective.zero_out(); - elec.pot->update_from_charge(elec.charge, &GlobalC::ucell); + elec.pot->update_from_charge(elec.charge, &ucell); v_effective = elec.pot->get_effective_v(); v_xc.create(PARAM.inp.nspin, rho_basis->nrxx); v_xc.zero_out(); const std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rho_basis->nrxx, elec.charge, &GlobalC::ucell); + = XC_Functional::v_xc(rho_basis->nrxx, elec.charge, &ucell); v_xc = std::get<2>(etxc_vtxc_v); GlobalC::paw_cell.calculate_force(v_effective.c, v_xc.c, rhor, force_paw); @@ -284,9 +285,9 @@ void Forces::cal_force(ModuleBase::matrix& force, double sum = 0.0; iat = 0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int ia = 0; ia < ucell.atoms[it].na; ia++) { force(iat, ipol) = forcelc(iat, ipol) + forceion(iat, ipol) + forcenl(iat, ipol) + forcecc(iat, ipol) + forcescc(iat, ipol); @@ -345,15 +346,15 @@ void Forces::cal_force(ModuleBase::matrix& force, ModuleBase::Mathzone::Cartesian_to_Direct(force(iat, 0), force(iat, 1), force(iat, 2), - GlobalC::ucell.a1.x, - GlobalC::ucell.a1.y, - GlobalC::ucell.a1.z, - GlobalC::ucell.a2.x, - GlobalC::ucell.a2.y, - GlobalC::ucell.a2.z, - GlobalC::ucell.a3.x, - GlobalC::ucell.a3.y, - GlobalC::ucell.a3.z, + ucell.a1.x, + ucell.a1.y, + ucell.a1.z, + ucell.a2.x, + ucell.a2.y, + ucell.a2.z, + ucell.a3.x, + ucell.a3.y, + ucell.a3.z, d1, d2, d3); @@ -368,15 +369,15 @@ void Forces::cal_force(ModuleBase::matrix& force, ModuleBase::Mathzone::Direct_to_Cartesian(force(iat, 0), force(iat, 1), force(iat, 2), - GlobalC::ucell.a1.x, - GlobalC::ucell.a1.y, - GlobalC::ucell.a1.z, - GlobalC::ucell.a2.x, - GlobalC::ucell.a2.y, - GlobalC::ucell.a2.z, - GlobalC::ucell.a3.x, - GlobalC::ucell.a3.y, - GlobalC::ucell.a3.z, + ucell.a1.x, + ucell.a1.y, + ucell.a1.z, + ucell.a2.x, + ucell.a2.y, + ucell.a2.z, + ucell.a3.x, + ucell.a3.y, + ucell.a3.z, d1, d2, d3); @@ -389,17 +390,17 @@ void Forces::cal_force(ModuleBase::matrix& force, GlobalV::ofs_running << std::setiosflags(std::ios::fixed) << std::setprecision(6) << std::endl; /*if(PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"LOCAL FORCE (Ry/Bohr)", forcelc); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"NONLOCAL FORCE (Ry/Bohr)", forcenl); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"NLCC FORCE (Ry/Bohr)", forcecc); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"ION FORCE (Ry/Bohr)", forceion); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"SCC FORCE (Ry/Bohr)", forcescc); - if(GlobalV::EFIELD) ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell,"EFIELD FORCE (Ry/Bohr)", + ModuleIO::print_force(GlobalV::ofs_running, ucell,"LOCAL FORCE (Ry/Bohr)", forcelc); + ModuleIO::print_force(GlobalV::ofs_running, ucell,"NONLOCAL FORCE (Ry/Bohr)", forcenl); + ModuleIO::print_force(GlobalV::ofs_running, ucell,"NLCC FORCE (Ry/Bohr)", forcecc); + ModuleIO::print_force(GlobalV::ofs_running, ucell,"ION FORCE (Ry/Bohr)", forceion); + ModuleIO::print_force(GlobalV::ofs_running, ucell,"SCC FORCE (Ry/Bohr)", forcescc); + if(GlobalV::EFIELD) ModuleIO::print_force(GlobalV::ofs_running, ucell,"EFIELD FORCE (Ry/Bohr)", force_e); }*/ /* - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell," TOTAL-FORCE (Ry/Bohr)", force); + ModuleIO::print_force(GlobalV::ofs_running, ucell," TOTAL-FORCE (Ry/Bohr)", force); if(INPUT.out_force) // pengfei 2016-12-20 { @@ -423,27 +424,27 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.test_force) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "LOCAL FORCE (eV/Angstrom)", forcelc, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "NONLOCAL FORCE (eV/Angstrom)", forcenl, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "NLCC FORCE (eV/Angstrom)", forcecc, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "ION FORCE (eV/Angstrom)", forceion, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "SCC FORCE (eV/Angstrom)", forcescc, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "LOCAL FORCE (eV/Angstrom)", forcelc, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "NONLOCAL FORCE (eV/Angstrom)", forcenl, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "NLCC FORCE (eV/Angstrom)", forcecc, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "ION FORCE (eV/Angstrom)", forceion, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "SCC FORCE (eV/Angstrom)", forcescc, false); if (PARAM.inp.use_paw) { ModuleIO::print_force(GlobalV::ofs_running, - GlobalC::ucell, + ucell, "PAW FORCE (eV/Angstrom)", forcepaw, false); } if (PARAM.inp.efield_flag) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD FORCE (eV/Angstrom)", force_e, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD FORCE (eV/Angstrom)", force_e, false); } if (PARAM.inp.gate_flag) { ModuleIO::print_force(GlobalV::ofs_running, - GlobalC::ucell, + ucell, "GATEFIELD FORCE (eV/Angstrom)", force_gate, false); @@ -451,19 +452,20 @@ void Forces::cal_force(ModuleBase::matrix& force, if (PARAM.inp.imp_sol) { ModuleIO::print_force(GlobalV::ofs_running, - GlobalC::ucell, + ucell, "IMP_SOL FORCE (eV/Angstrom)", forcesol, false); } } - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "TOTAL-FORCE (eV/Angstrom)", force, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "TOTAL-FORCE (eV/Angstrom)", force, false); ModuleBase::timer::tick("Forces", "cal_force"); return; } template -void Forces::cal_force_loc(ModuleBase::matrix& forcelc, +void Forces::cal_force_loc(const UnitCell& ucell, + ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& vloc, const Charge* const chr) @@ -515,11 +517,11 @@ void Forces::cal_force_loc(ModuleBase::matrix& forcelc, for (int iat = 0; iat < this->nat; ++iat) { // read `it` `ia` from the table - int it = GlobalC::ucell.iat2it[iat]; - int ia = GlobalC::ucell.iat2ia[iat]; + int it = ucell.iat2it[iat]; + int ia = ucell.iat2ia[iat]; for (int ig = 0; ig < rho_basis->npw; ig++) { - const double phase = ModuleBase::TWO_PI * (rho_basis->gcar[ig] * GlobalC::ucell.atoms[it].tau[ia]); + const double phase = ModuleBase::TWO_PI * (rho_basis->gcar[ig] * ucell.atoms[it].tau[ia]); double sinp, cosp; ModuleBase::libm::sincos(phase, &sinp, &cosp); const double factor @@ -528,9 +530,9 @@ void Forces::cal_force_loc(ModuleBase::matrix& forcelc, forcelc(iat, 1) += rho_basis->gcar[ig][1] * factor; forcelc(iat, 2) += rho_basis->gcar[ig][2] * factor; } - forcelc(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega); - forcelc(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega); - forcelc(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega); + forcelc(iat, 0) *= (ucell.tpiba * ucell.omega); + forcelc(iat, 1) *= (ucell.tpiba * ucell.omega); + forcelc(iat, 2) *= (ucell.tpiba * ucell.omega); } // this->print(GlobalV::ofs_running, "local forces", forcelc); @@ -541,7 +543,8 @@ void Forces::cal_force_loc(ModuleBase::matrix& forcelc, } template -void Forces::cal_force_ew(ModuleBase::matrix& forceion, +void Forces::cal_force_ew(const UnitCell& ucell, + ModuleBase::matrix& forceion, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf) { @@ -555,7 +558,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, blocking rho_basis->nrxnpwx for data locality. By blocking aux with block size 1024, - we can keep the blocked aux in L1 cache when iterating GlobalC::ucell.ntype loop + we can keep the blocked aux in L1 cache when iterating ucell.ntype loop performance will be better when number of atom is quite huge */ const int block_ig = 1024; @@ -571,7 +574,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, { aux[ig] = 0.0; } - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { double dzv; if (PARAM.inp.use_paw) @@ -582,7 +585,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, } else { - dzv = GlobalC::ucell.atoms[it].ncpp.zv; + dzv = ucell.atoms[it].ncpp.zv; } for (int ig = igb; ig < ig_end; ++ig) @@ -594,17 +597,17 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, // calculate total ionic charge double charge = 0.0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { if (PARAM.inp.use_paw) { #ifdef USE_PAW - charge += GlobalC::ucell.atoms[it].na * GlobalC::paw_cell.get_val(it); + charge += ucell.atoms[it].na * GlobalC::paw_cell.get_val(it); #endif } else { - charge += GlobalC::ucell.atoms[it].na * GlobalC::ucell.atoms[it].ncpp.zv; // mohan modify 2007-11-7 + charge += ucell.atoms[it].na * ucell.atoms[it].ncpp.zv; // mohan modify 2007-11-7 } } @@ -621,7 +624,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, ModuleBase::WARNING_QUIT("ewald", "Can't find optimal alpha."); } upperbound = 2.0 * charge * charge * sqrt(2.0 * alpha / ModuleBase::TWO_PI) - * erfc(sqrt(GlobalC::ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha)); + * erfc(sqrt(ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha)); } while (upperbound > 1.0e-6); #ifdef _OPENMP @@ -629,8 +632,8 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, #endif for (int ig = 0; ig < rho_basis->npw; ig++) { - aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * GlobalC::ucell.tpiba2 / alpha / 4.0) - / (rho_basis->gg[ig] * GlobalC::ucell.tpiba2); + aux[ig] *= ModuleBase::libm::exp(-1.0 * rho_basis->gg[ig] * ucell.tpiba2 / alpha / 4.0) + / (rho_basis->gg[ig] * ucell.tpiba2); } // set pos rho_basis->ig_gge0 to zero @@ -660,7 +663,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, int it_beg, ia_beg; ModuleBase::TASK_DIST_1D(num_threads, thread_id, this->nat, iat_beg, iat_end); iat_end = iat_beg + iat_end; - GlobalC::ucell.iat2iait(iat_beg, &ia_beg, &it_beg); + ucell.iat2iait(iat_beg, &ia_beg, &it_beg); int iat = iat_beg; int it = it_beg; @@ -686,9 +689,9 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, } else { - zv = GlobalC::ucell.atoms[it].ncpp.zv; + zv = ucell.atoms[it].ncpp.zv; } - it_fact = zv * ModuleBase::e2 * GlobalC::ucell.tpiba * ModuleBase::TWO_PI / GlobalC::ucell.omega * fact; + it_fact = zv * ModuleBase::e2 * ucell.tpiba * ModuleBase::TWO_PI / ucell.omega * fact; last_it = it; } @@ -696,7 +699,7 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, for (int ig = ig_beg; ig < ig_end; ig++) { const ModuleBase::Vector3 gcar = rho_basis->gcar[ig]; - const double arg = ModuleBase::TWO_PI * (gcar * GlobalC::ucell.atoms[it].tau[ia]); + const double arg = ModuleBase::TWO_PI * (gcar * ucell.atoms[it].tau[ia]); double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); double sumnb = -cosp * aux[ig].imag() + sinp * aux[ig].real(); @@ -715,13 +718,13 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, forceion(iat, 2) *= it_fact; ++iat; - GlobalC::ucell.step_iait(&ia, &it); + ucell.step_iait(&ia, &it); } // means that the processor contains G=0 term. if (rho_basis->ig_gge0 >= 0) { - double rmax = 5.0 / (sqrt(alpha) * GlobalC::ucell.lat0); + double rmax = 5.0 / (sqrt(alpha) * ucell.lat0); int nrm = 0; // output of rgen: the number of vectors in the sphere @@ -752,12 +755,12 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, if (iat1 != iat2) { ModuleBase::Vector3 d_tau - = GlobalC::ucell.atoms[T1].tau[I1] - GlobalC::ucell.atoms[T2].tau[I2]; - H_Ewald_pw::rgen(d_tau, rmax, irr, GlobalC::ucell.latvec, GlobalC::ucell.G, r, r2, nrm); + = ucell.atoms[T1].tau[I1] - ucell.atoms[T2].tau[I2]; + H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm); for (int n = 0; n < nrm; n++) { - const double rr = sqrt(r2[n]) * GlobalC::ucell.lat0; + const double rr = sqrt(r2[n]) * ucell.lat0; double factor; if (PARAM.inp.use_paw) @@ -766,15 +769,15 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, factor = GlobalC::paw_cell.get_val(T1) * GlobalC::paw_cell.get_val(T2) * ModuleBase::e2 / (rr * rr) * (erfc(sqa * rr) / rr + sq8a_2pi * ModuleBase::libm::exp(-alpha * rr * rr)) - * GlobalC::ucell.lat0; + * ucell.lat0; #endif } else { - factor = GlobalC::ucell.atoms[T1].ncpp.zv * GlobalC::ucell.atoms[T2].ncpp.zv + factor = ucell.atoms[T1].ncpp.zv * ucell.atoms[T2].ncpp.zv * ModuleBase::e2 / (rr * rr) * (erfc(sqa * rr) / rr + sq8a_2pi * ModuleBase::libm::exp(-alpha * rr * rr)) - * GlobalC::ucell.lat0; + * ucell.lat0; } forceion(iat1, 0) -= factor * r[n].x; @@ -783,10 +786,10 @@ void Forces::cal_force_ew(ModuleBase::matrix& forceion, } } ++iat2; - GlobalC::ucell.step_iait(&I2, &T2); + ucell.step_iait(&I2, &T2); } // atom b ++iat1; - GlobalC::ucell.step_iait(&I1, &T1); + ucell.step_iait(&I1, &T1); } // atom a delete[] r; diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h index eca52c95cf..c23f24f53e 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.h +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.h @@ -33,7 +33,8 @@ class Forces Forces(const int nat_in) : nat(nat_in){}; ~Forces(){}; - void cal_force(ModuleBase::matrix& force, + void cal_force(const UnitCell& ucell, + ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, @@ -48,11 +49,15 @@ class Forces int nat = 0; int npwx = 0; - void cal_force_loc(ModuleBase::matrix& forcelc, + void cal_force_loc(const UnitCell& ucell, + ModuleBase::matrix& forcelc, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& vloc, const Charge* const chr); - void cal_force_ew(ModuleBase::matrix& forceion, ModulePW::PW_Basis* rho_basis, const Structure_Factor* p_sf); + void cal_force_ew(const UnitCell& ucell, + ModuleBase::matrix& forceion, + ModulePW::PW_Basis* rho_basis, + const Structure_Factor* p_sf); void cal_force_cc(ModuleBase::matrix& forcecc, ModulePW::PW_Basis* rho_basis, const Charge* const chr, diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h b/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h index eff63d8c0b..461076eb4b 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_lcaopw.h @@ -16,15 +16,17 @@ namespace hamilt HamiltLIP(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, - pseudopot_cell_vnl* nlpp) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp){}; + pseudopot_cell_vnl* nlpp, + const UnitCell* ucell) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell){}; #ifdef __EXX HamiltLIP(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, + const UnitCell* ucell, Exx_Lip& exx_lip_in) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp), exx_lip(exx_lip_in){}; + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell), exx_lip(exx_lip_in){}; Exx_Lip& exx_lip; #endif }; diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp index abcc87ffc0..6272675398 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp @@ -22,15 +22,16 @@ template HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, - pseudopot_cell_vnl* nlpp) + pseudopot_cell_vnl* nlpp, + const UnitCell* ucell): ucell(ucell) { this->classname = "HamiltPW"; this->ppcell = nlpp; this->qq_nt = this->ppcell->template get_qq_nt_data(); this->qq_so = this->ppcell->template get_qq_so_data(); this->vkb = this->ppcell->template get_vkb_data(); - const auto tpiba2 = static_cast(GlobalC::ucell.tpiba2); - const auto tpiba = static_cast(GlobalC::ucell.tpiba); + const auto tpiba2 = static_cast(ucell->tpiba2); + const auto tpiba = static_cast(ucell->tpiba); const int* isk = pkv->isk.data(); const Real* gk2 = wfc_basis->get_gk2_data(); @@ -103,7 +104,7 @@ HamiltPW::HamiltPW(elecstate::Potential* pot_in, if (PARAM.inp.vnl_in_h) { Operator* nonlocal - = new Nonlocal>(isk, this->ppcell, &GlobalC::ucell, wfc_basis); + = new Nonlocal>(isk, this->ppcell, ucell, wfc_basis); if(this->ops == nullptr) { this->ops = nonlocal; @@ -290,9 +291,9 @@ void HamiltPW::sPsi(const T* psi_in, // psi // qq char transa = 'N'; char transb = 'N'; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell->ntype; it++) { - Atom* atoms = &GlobalC::ucell.atoms[it]; + Atom* atoms = &ucell->atoms[it]; if (atoms->ncpp.tvanp) { const int nh = atoms->ncpp.nh; @@ -309,7 +310,7 @@ void HamiltPW::sPsi(const T* psi_in, // psi } for (int ia = 0; ia < atoms->na; ia++) { - const int iat = GlobalC::ucell.itia2iat(it, ia); + const int iat = ucell->itia2iat(it, ia); gemm_op()(this->ctx, transa, transb, diff --git a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h index 7843c0d3f1..f87dca7745 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h +++ b/source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.h @@ -20,7 +20,7 @@ class HamiltPW : public Hamilt // otherwise return the real type of T(complex, complex) using Real = typename GetTypeReal::type; public: - HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp); + HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp,const UnitCell* ucell); template explicit HamiltPW(const HamiltPW* hamilt); ~HamiltPW(); @@ -38,6 +38,7 @@ class HamiltPW : public Hamilt protected: // used in sPhi, which are calculated in hPsi or sPhi const pseudopot_cell_vnl* ppcell = nullptr; + const UnitCell* const ucell = nullptr; mutable T* vkb = nullptr; Real* qq_nt = nullptr; T* qq_so = nullptr; diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/ekinetic_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/ekinetic_op.h index 7142fd4031..47024f47ba 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/ekinetic_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/ekinetic_op.h @@ -14,7 +14,7 @@ struct ekinetic_pw_op { /// \param nband : nbands /// \param npw : number of planewaves of current k point /// \param max_npw : max number of planewaves of all k points - /// \param tpiba2 : GlobalC::ucell.tpiba2 + /// \param tpiba2 : ucell.tpiba2 /// \param spin : current spin /// \param gk2_ik : GlobalC::wfcpw->gk2 /// \param tmpsi_in : intermediate array diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h index 539bdebe88..b9aaa6d468 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.h @@ -56,9 +56,9 @@ struct cal_force_nl_op /// @param forcenl_nc - the second dimension of matrix forcenl /// @param nbands - NBANDS /// @param nkb - number of k point - /// @param atom_nh - GlobalC::ucell.atoms[ii].ncpp.nh - /// @param atom_na - GlobalC::ucell.atoms[ii].na - /// @param tpiba - GlobalC::ucell.tpiba + /// @param atom_nh - ucell.atoms[ii].ncpp.nh + /// @param atom_na - ucell.atoms[ii].na + /// @param tpiba - ucell.tpiba /// @param d_wg - input parameter wg /// @param occ - if use the occupation of the bands /// @param d_ekb - input parameter ekb diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/meta_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/meta_op.h index d179add67e..293beeb8aa 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/meta_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/meta_op.h @@ -16,7 +16,7 @@ namespace hamilt { /// \param pol : loop of [0, 1, 2] /// \param npw : number of planewaves of current k point /// \param npwx : number of planewaves of all k points - /// \param tpiba : GlobalC::ucell.tpiba + /// \param tpiba : ucell.tpiba /// \param gcar : wfcpw->gcar /// \param kvec_c : wfcpw->kvec_c /// \param in : input hpsi diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h index 64cc01514a..af7d51523d 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.h @@ -23,7 +23,7 @@ struct cal_dbecp_noevc_nl_op /// @param npw - number of planewaves /// @param npwx - max number of planewaves /// @param ik - the current k point - /// @param tpiba - GlobalC::ucell.tpiba + /// @param tpiba - ucell.tpiba /// @param gcar - GlobalC::wfcpw->gcar /// @param kvec_c - GlobalC::wfcpw->kvec_c /// @param vkbi - _vkb0[ipol] @@ -69,8 +69,8 @@ struct cal_stress_nl_op /// @param deeq_2 - the second dimension of deeq /// @param deeq_3 - the third dimension of deeq /// @param deeq_4 - the forth dimension of deeq - /// @param atom_nh - GlobalC::ucell.atoms[ii].ncpp.nh - /// @param atom_na - GlobalC::ucell.atoms[ii].na + /// @param atom_nh - ucell.atoms[ii].ncpp.nh + /// @param atom_na - ucell.atoms[ii].na /// @param d_wg - input parameter wg /// @param occ - if use the occupation of the bands /// @param d_ekb - input parameter ekb diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.h index 3e02910abe..3f72cd1755 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.h @@ -20,11 +20,11 @@ struct cal_vnl_op /// @param npwx - number of planewaves of all k points /// @param tab_2 - the second dimension of the input table /// @param tab_3 - the third dimension of the input table - /// @param atom_nh - GlobalC::ucell.atoms[ii].ncpp.nh - /// @param atom_nb - GlobalC::ucell.atoms[it].ncpp.nbeta - /// @param atom_na - GlobalC::ucell.atoms[ii].na + /// @param atom_nh - ucell.atoms[ii].ncpp.nh + /// @param atom_nb - ucell.atoms[it].ncpp.nbeta + /// @param atom_na - ucell.atoms[ii].na /// @param DQ - PARAM.globalv.dq - /// @param tpiba - GlobalC::ucell.tpiba + /// @param tpiba - ucell.tpiba /// @param NEG_IMAG_UNIT - ModuleBase::NEG_IMAG_UNIT /// @param gk - GlobalC::wf.get_1qvec_cartesian /// @param ylm - the result of ModuleBase::YlmReal::Ylm_Real diff --git a/source/module_hamilt_pw/hamilt_pwdft/kernels/wf_op.h b/source/module_hamilt_pw/hamilt_pwdft/kernels/wf_op.h index 43ba561b3b..ae2a7ee105 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/kernels/wf_op.h +++ b/source/module_hamilt_pw/hamilt_pwdft/kernels/wf_op.h @@ -26,18 +26,18 @@ struct cal_sk_op { /// @param eigts1_nc - GlobalC::sf.eigts1.nc /// @param eigts2_nc - GlobalC::sf.eigts2.nc /// @param eigts3_nc - GlobalC::sf.eigts3.nc - /// @param atom_na - GlobalC::ucell.atoms[ii].na + /// @param atom_na - ucell.atoms[ii].na /// @param igl2isz - wfc_basis->igl2isz_k /// @param is2fftixy - wfc_basis->is2fftixy /// @param TWO_PI - ModuleBase::TWO_PI /// @param kvec_c - GlobalC::kv.kvec_c - /// @param atom_tau - GlobalC::ucell.atoms[it].tau + /// @param atom_tau - ucell.atoms[it].tau /// @param eigts1 - GlobalC::sf.eigts1 /// @param eigts2 - GlobalC::sf.eigts2 /// @param eigts3 - GlobalC::sf.eigts3 /// /// Output Parameters - /// @param sk - output results matrix with size GlobalC::ucell.nat * npw + /// @param sk - output results matrix with size ucell.nat * npw void operator()(const Device* ctx, const int& ik, const int& ntype, diff --git a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/nonlocal_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/nonlocal_pw.cpp index 96e3db88ee..b8d9a3f086 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/nonlocal_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/nonlocal_pw.cpp @@ -46,7 +46,7 @@ void Nonlocal>::init(const int ik_in) // Calculate nonlocal pseudopotential vkb if(this->ppcell->nkb > 0) //xiaohui add 2013-09-02. Attention... { - this->ppcell->getvnl(this->ctx, this->ik, this->vkb); + this->ppcell->getvnl(this->ctx, *this->ucell, this->ik, this->vkb); } if(this->next_op != nullptr) diff --git a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.cpp index efe4493831..858e6b3fd5 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.cpp @@ -32,7 +32,7 @@ void Velocity::init(const int ik_in) // Calculate nonlocal pseudopotential vkb if(this->ppcell->nkb > 0 && this->nonlocal) { - this->ppcell->getgradq_vnl(ik_in); + this->ppcell->getgradq_vnl(*this->ucell,ik_in); } } diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h index fbd51930e3..6d5ee5581e 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func.h +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func.h @@ -70,19 +70,22 @@ class Stress_Func const psi::Psi, Device>* psi_in = nullptr); // electron kinetic part in PW basis // 2) the stress from the Hartree term - void stress_har(ModuleBase::matrix& sigma, + void stress_har(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw, const Charge* const chr); // hartree part in PW or LCAO basis // 3) the stress from the ewald term (ion-ion intraction under // periodic boundary conditions). - void stress_ewa(ModuleBase::matrix& sigma, + void stress_ewa(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw); // ewald part in PW or LCAO basis // 4) the stress from the local pseudopotentials - void stress_loc(ModuleBase::matrix& sigma, + void stress_loc(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& vloc, const Structure_Factor* p_sf, @@ -103,7 +106,8 @@ class Stress_Func * D V(g^2) / D g^2 = 4pi e^2/omegai /G^4 * */ - void dvloc_coulomb(const FPTYPE& zp, + void dvloc_coulomb(const UnitCell& ucell, + const FPTYPE& zp, FPTYPE* dvloc, ModulePW::PW_Basis* rho_basis); // used in local pseudopotential stress @@ -125,10 +129,12 @@ class Stress_Func int type); // used in nonlinear core correction stress // 6) the stress from the exchange-correlation functional term - void stress_gga(ModuleBase::matrix& sigma, + void stress_gga(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Charge* const chr); // gga part in both PW and LCAO basis - void stress_mgga(ModuleBase::matrix& sigma, + void stress_mgga(const UnitCell& ucell, + ModuleBase::matrix& sigma, const ModuleBase::matrix& wg, const ModuleBase::matrix& v_ofk, const Charge* const chr, diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp index d948a8c19b..fd6bf8dc0a 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp @@ -11,15 +11,18 @@ //calcualte the Ewald stress term in PW and LCAO template -void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw) +void Stress_Func::stress_ewa(const UnitCell& ucell, + ModuleBase::matrix& sigma, + ModulePW::PW_Basis* rho_basis, + const bool is_pw) { ModuleBase::TITLE("Stress_Func","stress_ewa"); ModuleBase::timer::tick("Stress_Func","stress_ewa"); FPTYPE charge=0; - for(int it=0; it < GlobalC::ucell.ntype; it++) + for(int it=0; it < ucell.ntype; it++) { - charge = charge + GlobalC::ucell.atoms[it].ncpp.zv * GlobalC::ucell.atoms[it].na; + charge = charge + ucell.atoms[it].ncpp.zv * ucell.atoms[it].na; } //choose alpha in order to have convergence in the sum over G //upperbound is a safe upper bound for the error ON THE ENERGY @@ -30,7 +33,7 @@ void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW alpha-=0.1; if(alpha==0.0) ModuleBase::WARNING_QUIT("stres_ew", "optimal alpha not found"); - upperbound =ModuleBase::e2 * pow(charge,2) * sqrt( 2 * alpha / (ModuleBase::TWO_PI)) * erfc(sqrt(GlobalC::ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha)); + upperbound =ModuleBase::e2 * pow(charge,2) * sqrt( 2 * alpha / (ModuleBase::TWO_PI)) * erfc(sqrt(ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha)); } while(upperbound>1e-7); @@ -40,7 +43,7 @@ void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW const int ig0 = rho_basis->ig_gge0; if( ig0 >= 0) { - sdewald = (ModuleBase::TWO_PI) * ModuleBase::e2 / 4.0 / alpha * pow(charge/GlobalC::ucell.omega,2); + sdewald = (ModuleBase::TWO_PI) * ModuleBase::e2 / 4.0 / alpha * pow(charge/ucell.omega,2); } else { @@ -79,27 +82,27 @@ void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW for(; ig < ig_end; ig++) { if(ig == ig0) continue; - g2 = rho_basis->gg[ig]* GlobalC::ucell.tpiba2; + g2 = rho_basis->gg[ig]* ucell.tpiba2; g2a = g2 /4.0/alpha; rhostar=std::complex(0.0,0.0); - for(int it=0; it < GlobalC::ucell.ntype; it++) + for(int it=0; it < ucell.ntype; it++) { - for(int i=0; igcar[ig] * GlobalC::ucell.atoms[it].tau[i]) * (ModuleBase::TWO_PI); + arg = (rho_basis->gcar[ig] * ucell.atoms[it].tau[i]) * (ModuleBase::TWO_PI); FPTYPE sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); - rhostar = rhostar + std::complex(GlobalC::ucell.atoms[it].ncpp.zv * cosp,GlobalC::ucell.atoms[it].ncpp.zv * sinp); + rhostar = rhostar + std::complex(ucell.atoms[it].ncpp.zv * cosp,ucell.atoms[it].ncpp.zv * sinp); } } - rhostar /= GlobalC::ucell.omega; + rhostar /= ucell.omega; sewald = fact* (ModuleBase::TWO_PI) * ModuleBase::e2 * ModuleBase::libm::exp(-g2a) / g2 * pow(std::abs(rhostar),2); local_sdewald -= sewald; for(int l=0;l<3;l++) { for(int m=0;mgcar[ig][l] * rho_basis->gcar[ig][m] / g2 * (g2a + 1); + local_sigma(l, m) += sewald * ucell.tpiba2 * 2.0 * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / g2 * (g2a + 1); } } } @@ -123,25 +126,25 @@ void Stress_Func::stress_ewa(ModuleBase::matrix& sigma, ModulePW FPTYPE sqa = sqrt(alpha); FPTYPE sq8a_2pi = sqrt(8 * alpha / (ModuleBase::TWO_PI)); - rmax = 4.0/sqa/GlobalC::ucell.lat0; + rmax = 4.0/sqa/ucell.lat0; // collapse it, ia, jt, ja loop into a single loop long long ijat, ijat_end; int it, i, jt, j; - ModuleBase::TASK_DIST_1D(num_threads, thread_id, (long long)GlobalC::ucell.nat * GlobalC::ucell.nat, ijat, ijat_end); + ModuleBase::TASK_DIST_1D(num_threads, thread_id, (long long)ucell.nat * ucell.nat, ijat, ijat_end); ijat_end = ijat + ijat_end; - GlobalC::ucell.ijat2iaitjajt(ijat, &i, &it, &j, &jt); + ucell.ijat2iaitjajt(ijat, &i, &it, &j, &jt); while (ijat < ijat_end) { //calculate tau[na]-tau[nb] - d_tau = GlobalC::ucell.atoms[it].tau[i] - GlobalC::ucell.atoms[jt].tau[j]; + d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j]; //generates nearest-neighbors shells - H_Ewald_pw::rgen(d_tau, rmax, irr, GlobalC::ucell.latvec, GlobalC::ucell.G, r, r2, nrm); + H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm); for(int nr=0 ; nr::stress_ewa(ModuleBase::matrix& sigma, ModulePW }//end nr ++ijat; - GlobalC::ucell.step_jajtiait(&j, &jt, &i, &it); + ucell.step_jajtiait(&j, &jt, &i, &it); } delete[] r; diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_gga.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_gga.cpp index e5b261fd15..54c2ee3b51 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_gga.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_gga.cpp @@ -5,7 +5,8 @@ //calculate the GGA stress correction in PW and LCAO template -void Stress_Func::stress_gga(ModuleBase::matrix& sigma, +void Stress_Func::stress_gga(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const Charge* const chr) { @@ -26,7 +27,7 @@ void Stress_Func::stress_gga(ModuleBase::matrix& sigma, // call gradcorr to evaluate gradient correction to stress // the first three terms are etxc, vtxc and v, which // is not used here, so dummy variables are used. - XC_Functional::gradcorr(dum1, dum2, dum3, chr, rho_basis, &GlobalC::ucell, stress_gga, 1); + XC_Functional::gradcorr(dum1, dum2, dum3, chr, rho_basis, &ucell, stress_gga, 1); for(int l = 0;l< 3;l++) { diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_har.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_har.cpp index 9858b0186e..9f2e702147 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_har.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_har.cpp @@ -6,7 +6,11 @@ //calculate the Hartree part in PW or LCAO base template -void Stress_Func::stress_har(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw, const Charge* const chr) +void Stress_Func::stress_har(const UnitCell& ucell, + ModuleBase::matrix& sigma, + ModulePW::PW_Basis* rho_basis, + const bool is_pw, + const Charge* const chr) { ModuleBase::TITLE("Stress_Func","stress_har"); ModuleBase::timer::tick("Stress_Func","stress_har"); @@ -65,10 +69,10 @@ void Stress_Func::stress_har(ModuleBase::matrix& sigma, ModulePW const FPTYPE g2 = rho_basis->gg[ig]; if(g2 < 1e-8) { continue; } - //const FPTYPE fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (GlobalC::ucell.tpiba2 * GlobalC::sf.gg [ig]); + //const FPTYPE fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (ucell.tpiba2 * GlobalC::sf.gg [ig]); //ehart += ( conj( Porter[j] ) * Porter[j] ).real() * fac; //vh_g[ig] = fac * Porter[j]; - FPTYPE shart= ( conj( aux[ig] ) * aux[ig] ).real()/(GlobalC::ucell.tpiba2 * g2); + FPTYPE shart= ( conj( aux[ig] ) * aux[ig] ).real()/(ucell.tpiba2 * g2); for(int l=0;l<3;l++) { for(int m=0;m::stress_har(ModuleBase::matrix& sigma, ModulePW } // Parallel_Reduce::reduce_pool( ehart ); -// ehart *= 0.5 * GlobalC::ucell.omega; +// ehart *= 0.5 * ucell.omega; //psic(:)=(0.0,0.0) if(is_pw&&PARAM.globalv.gamma_only_pw) { @@ -125,8 +129,8 @@ void Stress_Func::stress_har(ModuleBase::matrix& sigma, ModulePW for(int l=0;l<3;l++) { - if(is_pw) { sigma(l,l) -= elecstate::H_Hartree_pw::hartree_energy /GlobalC::ucell.omega; - } else { sigma(l,l) += elecstate::H_Hartree_pw::hartree_energy /GlobalC::ucell.omega; + if(is_pw) { sigma(l,l) -= elecstate::H_Hartree_pw::hartree_energy /ucell.omega; + } else { sigma(l,l) += elecstate::H_Hartree_pw::hartree_energy /ucell.omega; } for(int m=0;m -void Stress_Func::stress_loc(ModuleBase::matrix& sigma, +void Stress_Func::stress_loc(const UnitCell& ucell, + ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const ModuleBase::matrix& vloc, const Structure_Factor* p_sf, @@ -69,7 +70,7 @@ void Stress_Func::stress_loc(ModuleBase::matrix& sigma, if(is_pw) { #pragma omp parallel for collapse(2) reduction(+:evloc) - for (int it=0; itnpw; ig++) { @@ -83,15 +84,15 @@ void Stress_Func::stress_loc(ModuleBase::matrix& sigma, } } } - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < ucell.ntype; ++it) { - const Atom* atom = &GlobalC::ucell.atoms[it]; + const Atom* atom = &ucell.atoms[it]; if(atom->coulomb_potential) { // // special case: pseudopotential is coulomb 1/r potential // - this->dvloc_coulomb (atom->ncpp.zv, dvloc.data(), rho_basis); + this->dvloc_coulomb (ucell,atom->ncpp.zv, dvloc.data(), rho_basis); // } else @@ -100,7 +101,7 @@ void Stress_Func::stress_loc(ModuleBase::matrix& sigma, // normal case: dvloc contains dV_loc(G)/dG // this->dvloc_of_g ( atom->ncpp.msh, atom->ncpp.rab.data(), atom->ncpp.r.data(), - atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, GlobalC::ucell); + atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis, ucell); // } #ifndef _OPENMP @@ -119,7 +120,7 @@ void Stress_Func::stress_loc(ModuleBase::matrix& sigma, { local_sigma(l, m) = local_sigma(l, m) + (conj(aux[ig]) * p_sf->strucFac(it, ig)).real() * 2.0 - * dvloc[rho_basis->ig2igg[ig]] * GlobalC::ucell.tpiba2 + * dvloc[rho_basis->ig2igg[ig]] * ucell.tpiba2 * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] * fact; } } @@ -277,7 +278,10 @@ const UnitCell& ucell_in } template -void Stress_Func::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc, ModulePW::PW_Basis* rho_basis) +void Stress_Func::dvloc_coulomb(const UnitCell& ucell, + const FPTYPE& zp, + FPTYPE* dvloc, + ModulePW::PW_Basis* rho_basis) { int igl0; // start from |G|=0 or not. @@ -295,8 +299,8 @@ void Stress_Func::dvloc_coulomb(const FPTYPE& zp, FPTYPE* dvloc, #endif for (int i = igl0; i < rho_basis->ngg; i++) { - dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / GlobalC::ucell.omega - / pow((GlobalC::ucell.tpiba2 * rho_basis->gg_uniq[i]), 2); + dvloc[i] = ModuleBase::FOUR_PI * zp * ModuleBase::e2 / ucell.omega + / pow((ucell.tpiba2 * rho_basis->gg_uniq[i]), 2); } return; diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp index 352634bee1..85f49d1138 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_mgga.cpp @@ -9,7 +9,8 @@ // calculate the Pulay term of mGGA stress correction in PW template -void Stress_Func::stress_mgga(ModuleBase::matrix& sigma, +void Stress_Func::stress_mgga(const UnitCell& ucell, + ModuleBase::matrix& sigma, const ModuleBase::matrix& wg, const ModuleBase::matrix& v_ofk, const Charge* const chr, @@ -52,9 +53,9 @@ void Stress_Func::stress_mgga(ModuleBase::matrix& sigma, for (int ibnd = 0; ibnd < PARAM.inp.nbands; ibnd++) { - const FPTYPE w1 = wg(ik, ibnd) / GlobalC::ucell.omega; + const FPTYPE w1 = wg(ik, ibnd) / ucell.omega; const std::complex* psi = &psi_in[0](ik, ibnd, 0); - XC_Functional::grad_wfc, Device>(ik, GlobalC::ucell.tpiba, wfc_basis, psi, gradwfc.data>()); + XC_Functional::grad_wfc, Device>(ik, ucell.tpiba, wfc_basis, psi, gradwfc.data>()); cal_stress_mgga_solver( current_spin, nrxx, w1, gradwfc.data>(), crosstaus.data()); } // band loop diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp index 713a233990..2a6925c912 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp @@ -66,20 +66,21 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->stress_kin(sigmakin, this->pelec->wg, p_symm, p_kv, wfc_basis, ucell, d_psi_in); // hartree contribution - this->stress_har(sigmahar, rho_basis, 1, pelec->charge); + this->stress_har(ucell,sigmahar, rho_basis, 1, pelec->charge); // ewald contribution - this->stress_ewa(sigmaewa, rho_basis, 1); + this->stress_ewa(ucell,sigmaewa, rho_basis, 1); // xc contribution: add gradient corrections(non diagonal) for (int i = 0; i < 3; i++) { sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega; } - this->stress_gga(sigmaxc, rho_basis, pelec->charge); + this->stress_gga(ucell,sigmaxc, rho_basis, pelec->charge); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - this->stress_mgga(sigmaxc, + this->stress_mgga(ucell, + sigmaxc, this->pelec->wg, this->pelec->pot->get_effective_vofk(), pelec->charge, @@ -89,7 +90,7 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, } // local contribution - this->stress_loc(sigmaloc, rho_basis, nlpp.vloc, p_sf, 1, pelec->charge); + 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); diff --git a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp index 6dd2d845f9..0fad7014bd 100644 --- a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.cpp @@ -10,10 +10,11 @@ HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, + const UnitCell* ucell, const int& npol, double* emin_in, double* emax_in) - : HamiltPW(pot_in, wfc_basis, p_kv, nlpp), ngk(p_kv->ngk) + : HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell), ngk(p_kv->ngk) { this->classname = "HamiltSdftPW"; this->npwk_max = wfc_basis->npwk_max; diff --git a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h index 599060dca4..709c87cfa4 100644 --- a/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h +++ b/source/module_hamilt_pw/hamilt_stodft/hamilt_sdft_pw.h @@ -24,6 +24,7 @@ class HamiltSdftPW : public HamiltPW ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, + const UnitCell* ucell, const int& npol, double* emin_in, double* emax_in); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp index e3289273a6..93ad9128e0 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp @@ -41,8 +41,8 @@ void Sto_Forces::cal_stoforce(ModuleBase::matrix& force, ModuleBase::matrix forcecc(this->nat, 3); ModuleBase::matrix forcenl(this->nat, 3); ModuleBase::matrix forcescc(this->nat, 3); - this->cal_force_loc(forcelc, rho_basis, nlpp.vloc, chr); - this->cal_force_ew(forceion, rho_basis, p_sf); + this->cal_force_loc(ucell,forcelc, rho_basis, nlpp.vloc, chr); + this->cal_force_ew(ucell,forceion, rho_basis, p_sf); this->cal_sto_force_nl(forcenl, wg, pkv, wfc_basis, p_sf, nlpp, ucell, psi, stowf); this->cal_force_cc(forcecc, rho_basis, chr, nlpp.numeric, ucell); this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, nlpp.numeric, ucell); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp index 08f2b894e1..bcfbd2da61 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_iter.cpp @@ -560,7 +560,8 @@ void Stochastic_Iter::sum_stoeband(Stochastic_WF& stowf, } template -void Stochastic_Iter::cal_storho(Stochastic_WF& stowf, +void Stochastic_Iter::cal_storho(const UnitCell& ucell, + Stochastic_WF& stowf, elecstate::ElecStatePW* pes, ModulePW::PW_Basis_K* wfc_basis) { @@ -651,12 +652,12 @@ void Stochastic_Iter::cal_storho(Stochastic_WF& stowf, #endif for (int ir = 0; ir < nrxx; ++ir) { - sto_rho[is][ir] /= GlobalC::ucell.omega; + sto_rho[is][ir] /= ucell.omega; sto_ne += sto_rho[is][ir]; } } - sto_ne *= GlobalC::ucell.omega / wfc_basis->nxyz; + sto_ne *= ucell.omega / wfc_basis->nxyz; #ifdef __MPI MPI_Allreduce(MPI_IN_PLACE, &sto_ne, 1, MPI_DOUBLE, MPI_SUM, POOL_WORLD); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_iter.h b/source/module_hamilt_pw/hamilt_stodft/sto_iter.h index f990b0078b..76cff4a0d9 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_iter.h +++ b/source/module_hamilt_pw/hamilt_stodft/sto_iter.h @@ -60,11 +60,13 @@ class Stochastic_Iter /** * @brief calculate the density * + * @param ucell reference to unit cell * @param stowf stochastic wave function * @param pes elecstate * @param wfc_basis wfc pw basis */ - void cal_storho(Stochastic_WF& stowf, + void cal_storho(const UnitCell& ucell, + Stochastic_WF& stowf, elecstate::ElecStatePW* pes, ModulePW::PW_Basis_K* wfc_basis); 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 0ca1d4dce2..f44172c52a 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp @@ -39,20 +39,20 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->sto_stress_kin(sigmakin, wg, p_symm, p_kv, wfc_basis, psi_in, stowf); // hartree contribution - this->stress_har(sigmahar, rho_basis, true, chr); + this->stress_har(ucell_in,sigmahar, rho_basis, true, chr); // ewald contribution - this->stress_ewa(sigmaewa, rho_basis, true); + this->stress_ewa(ucell_in,sigmaewa, rho_basis, true); // xc contribution: add gradient corrections(non diagonal) for (int i = 0; i < 3; ++i) { sigmaxc(i, i) = -(elec.f_en.etxc - elec.f_en.vtxc) / this->ucell->omega; } - this->stress_gga(sigmaxc, rho_basis, chr); + this->stress_gga(ucell_in,sigmaxc, rho_basis, chr); // local contribution - this->stress_loc(sigmaloc, rho_basis, nlpp->vloc, p_sf, true, chr); + 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); diff --git a/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp b/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp index bdbad81516..305cc548f9 100644 --- a/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/test/test_hamilt_sto.cpp @@ -11,7 +11,7 @@ void elecstate::Potential::cal_v_eff(Charge const*, UnitCell const*, ModuleBase: void elecstate::Potential::cal_fixed_v(double*){} template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl*){} +hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl*,const UnitCell*){} template hamilt::HamiltPW::~HamiltPW(){ delete this->ops; @@ -70,7 +70,7 @@ class TestHamiltSto : public ::testing::Test p_kv = new K_Vectors(); std::vector ngk = {2}; p_kv->ngk = ngk; - hamilt_sto = new hamilt::HamiltSdftPW, base_device::DEVICE_CPU>(pot, wfc_basis, p_kv, nullptr, npol, &emin, &emax); + hamilt_sto = new hamilt::HamiltSdftPW, base_device::DEVICE_CPU>(pot, wfc_basis, p_kv, nullptr, nullptr, npol, &emin, &emax); hamilt_sto->ops = new TestOp, base_device::DEVICE_CPU>(); } diff --git a/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp index d947801700..2e429dfa3f 100644 --- a/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp @@ -8,7 +8,7 @@ ***********************************************/ template -hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv,pseudopot_cell_vnl*){} +hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv,pseudopot_cell_vnl*,const UnitCell*){} template hamilt::HamiltPW::~HamiltPW(){}; template @@ -21,9 +21,10 @@ hamilt::HamiltSdftPW::HamiltSdftPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* p_kv, pseudopot_cell_vnl* nlpp, + const UnitCell* ucell, const int& npol, double* emin_in, - double* emax_in): HamiltPW(pot_in, wfc_basis, p_kv, nlpp), ngk(p_kv->ngk){} + double* emax_in): HamiltPW(pot_in, wfc_basis, p_kv, nlpp,ucell), ngk(p_kv->ngk){} template void hamilt::HamiltSdftPW::hPsi_norm(const T* psi_in, T* hpsi, const int& nbands){} diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp index 2eec4fa6ca..670b723b78 100644 --- a/source/module_hsolver/hsolver_pw_sdft.cpp +++ b/source/module_hsolver/hsolver_pw_sdft.cpp @@ -11,7 +11,8 @@ namespace hsolver { template -void HSolverPW_SDFT::solve(hamilt::Hamilt* pHamilt, +void HSolverPW_SDFT::solve(const UnitCell& ucell, + hamilt::Hamilt* pHamilt, psi::Psi& psi, psi::Psi& psi_cpu, elecstate::ElecState* pes, @@ -141,7 +142,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt* pHamilt, pes_pw->psiToRho(psi); } // calculate stochastic rho - stoiter.cal_storho(stowf, pes_pw, wfc_basis); + stoiter.cal_storho(ucell, stowf, pes_pw,wfc_basis); // will do rho symmetry and energy calculation in esolver ModuleBase::timer::tick("HSolverPW_SDFT", "solve"); diff --git a/source/module_hsolver/hsolver_pw_sdft.h b/source/module_hsolver/hsolver_pw_sdft.h index 0089c16e23..e47b77afe4 100644 --- a/source/module_hsolver/hsolver_pw_sdft.h +++ b/source/module_hsolver/hsolver_pw_sdft.h @@ -42,7 +42,8 @@ class HSolverPW_SDFT : public HSolverPW stoiter.init(pkv, wfc_basis_in, stowf, stoche, p_hamilt_sto); } - void solve(hamilt::Hamilt* pHamilt, + void solve(const UnitCell& ucell, + hamilt::Hamilt* pHamilt, psi::Psi& psi, psi::Psi& psi_cpu, elecstate::ElecState* pes, diff --git a/source/module_hsolver/test/diago_bpcg_test.cpp b/source/module_hsolver/test/diago_bpcg_test.cpp index 6274f75d74..d15b539c7d 100644 --- a/source/module_hsolver/test/diago_bpcg_test.cpp +++ b/source/module_hsolver/test/diago_bpcg_test.cpp @@ -97,7 +97,7 @@ class DiagoBPCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_cg_float_test.cpp b/source/module_hsolver/test/diago_cg_float_test.cpp index bab85fb179..1e33d67e14 100644 --- a/source/module_hsolver/test/diago_cg_float_test.cpp +++ b/source/module_hsolver/test/diago_cg_float_test.cpp @@ -108,7 +108,7 @@ class DiagoCGPrepare float *en = new float[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_cg_test.cpp b/source/module_hsolver/test/diago_cg_test.cpp index ada6b90210..f9b138efd9 100644 --- a/source/module_hsolver/test/diago_cg_test.cpp +++ b/source/module_hsolver/test/diago_cg_test.cpp @@ -104,7 +104,7 @@ class DiagoCGPrepare double *en = new double[npw]; int ik = 1; hamilt::Hamilt>* ha; - ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); + ha =new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); int* ngk = new int [1]; //psi::Psi> psi(ngk,ik,nband,npw); psi::Psi> psi; diff --git a/source/module_hsolver/test/diago_david_float_test.cpp b/source/module_hsolver/test/diago_david_float_test.cpp index 98a5d7c775..f6d97f87c9 100644 --- a/source/module_hsolver/test/diago_david_float_test.cpp +++ b/source/module_hsolver/test/diago_david_float_test.cpp @@ -82,7 +82,7 @@ class DiagoDavPrepare //do Diago_David::diag() float* en = new float[npw]; hamilt::Hamilt> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/module_hsolver/test/diago_david_real_test.cpp b/source/module_hsolver/test/diago_david_real_test.cpp index 67ca3fde8a..f6597fa78f 100644 --- a/source/module_hsolver/test/diago_david_real_test.cpp +++ b/source/module_hsolver/test/diago_david_real_test.cpp @@ -81,7 +81,7 @@ class DiagoDavPrepare //do Diago_David::diag() double* en = new double[npw]; hamilt::Hamilt* phm; - phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW(nullptr, nullptr, nullptr, nullptr,nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; diff --git a/source/module_hsolver/test/diago_david_test.cpp b/source/module_hsolver/test/diago_david_test.cpp index 237233a941..de88aadfe0 100644 --- a/source/module_hsolver/test/diago_david_test.cpp +++ b/source/module_hsolver/test/diago_david_test.cpp @@ -48,8 +48,10 @@ void lapackEigen(int& npw, std::vector>& hm, double* e, boo char tmp_c1 = 'V', tmp_c2 = 'U'; zheev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, rwork, &info); end = clock(); - if(info) std::cout << "ERROR: Lapack solver, info=" << info <> *phm; - phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr); + phm = new hamilt::HamiltPW>(nullptr, nullptr, nullptr, nullptr,nullptr); #ifdef __MPI const hsolver::diag_comm_info comm_info = {MPI_COMM_WORLD, mypnum, nprocs}; @@ -113,7 +116,7 @@ class DiagoDavPrepare const int ld_psi, const int nvec) { auto psi_iter_wrapper = psi::Psi>(psi_in, 1, nvec, ld_psi, nullptr); - psi::Range bands_range(1, 0, 0, nvec-1); + psi::Range bands_range(true, 0, 0, nvec-1); using hpsi_info = typename hamilt::Operator>::hpsi_info; hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out); phm->ops->hPsi(info); @@ -133,7 +136,8 @@ class DiagoDavPrepare if(mypnum == 0) { - if (DETAILINFO) std::cout<<"diag Run time: "<< use_time << std::endl; + if (DETAILINFO) { std::cout<<"diag Run time: "<< use_time << std::endl; +} for(int i=0;i {}; TEST_P(DiagoDavTest,RandomHamilt) { DiagoDavPrepare ddp = GetParam(); - if (DETAILINFO&&ddp.mypnum==0) std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" + if (DETAILINFO&&ddp.mypnum==0) { std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity=" << ddp.sparsity << ", eps=" << ddp.eps << std::endl; +} HPsi> hpsi(ddp.nband, ddp.npw, ddp.sparsity); DIAGOTEST::hmatrix = hpsi.hamilt(); @@ -238,7 +243,8 @@ int main(int argc, char **argv) testing::InitGoogleTest(&argc, argv); ::testing::TestEventListeners &listeners = ::testing::UnitTest::GetInstance()->listeners(); - if (myrank != 0) delete listeners.Release(listeners.default_result_printer()); + if (myrank != 0) { delete listeners.Release(listeners.default_result_printer()); +} int result = RUN_ALL_TESTS(); if (myrank == 0 && result != 0) diff --git a/source/module_hsolver/test/diago_mock.h b/source/module_hsolver/test/diago_mock.h index 973d01593b..e63022f43d 100644 --- a/source/module_hsolver/test/diago_mock.h +++ b/source/module_hsolver/test/diago_mock.h @@ -574,7 +574,7 @@ template<> void hamilt::HamiltPW::updateHk(const int ik) return; } -template<> hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) +template<> hamilt::HamiltPW::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) { this->ops = new OperatorMock; } @@ -589,7 +589,7 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) +template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) { this->ops = new OperatorMock>; } @@ -604,7 +604,7 @@ template<> void hamilt::HamiltPW>::updateHk(const int ik) return; } -template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*) +template<> hamilt::HamiltPW>::HamiltPW(elecstate::Potential* pot_in, ModulePW::PW_Basis_K* wfc_basis, K_Vectors* pkv, pseudopot_cell_vnl*,const UnitCell*) { this->ops = new OperatorMock>; } diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp index 81e09579cb..4d3db33096 100644 --- a/source/module_hsolver/test/test_hsolver_sdft.cpp +++ b/source/module_hsolver/test/test_hsolver_sdft.cpp @@ -152,7 +152,8 @@ void Stochastic_Iter::sum_stoeband(Stochastic_WF& stowf, } template -void Stochastic_Iter::cal_storho(Stochastic_WF& stowf, +void Stochastic_Iter::cal_storho(const UnitCell& ucell, + Stochastic_WF& stowf, elecstate::ElecStatePW* pes, ModulePW::PW_Basis_K* wfc_basis) {