Skip to content

Commit

Permalink
Fix: GlobalC in Sto_stress_PW
Browse files Browse the repository at this point in the history
  • Loading branch information
dyzheng committed May 13, 2024
1 parent 1db1b36 commit 9e9d748
Show file tree
Hide file tree
Showing 4 changed files with 35 additions and 24 deletions.
4 changes: 3 additions & 1 deletion source/module_esolver/esolver_sdft_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,9 @@ void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress)
pw_wfc,
this->psi,
this->stowf,
pelec->charge);
pelec->charge,
&GlobalC::ppcell,
GlobalC::ucell);
}


Expand Down
1 change: 1 addition & 0 deletions source/module_hamilt_pw/hamilt_pwdft/stress_func.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -299,6 +299,7 @@ class Stress_Func
using cal_vkb_op = hamilt::cal_vkb_op<FPTYPE, Device>;
using cal_vkb_deri_op = hamilt::cal_vkb_deri_op<FPTYPE, Device>;

protected:
pseudopot_cell_vnl* nlpp = nullptr;
const UnitCell* ucell = nullptr;

Expand Down
47 changes: 26 additions & 21 deletions source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,14 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
ModulePW::PW_Basis_K* wfc_basis,
const psi::Psi<complex<double>>* psi_in,
Stochastic_WF& stowf,
const Charge* const chr)
const Charge* const chr,
pseudopot_cell_vnl* nlpp_in,
const UnitCell& ucell_in)
{
ModuleBase::TITLE("Sto_Stress_PW", "cal_stress");
ModuleBase::timer::tick("Sto_Stress_PW", "cal_stress");
const ModuleBase::matrix& wg = elec.wg;
this->ucell = &ucell_in;
sigmatot.create(3, 3);
ModuleBase::matrix sigmaxc(3, 3);
ModuleBase::matrix sigmahar(3, 3);
Expand All @@ -40,7 +43,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
// 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) / GlobalC::ucell.omega;
sigmaxc(i, i) = -(elec.f_en.etxc - elec.f_en.vtxc) / this->ucell->omega;
}
stress_gga(sigmaxc, rho_basis, chr);

Expand All @@ -51,7 +54,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
stress_cc(sigmaxcc, rho_basis, p_sf, 1, chr);

// nonlocal
sto_stress_nl(sigmanl, wg, p_sf, p_symm, p_kv, wfc_basis, psi_in, stowf);
sto_stress_nl(sigmanl, wg, p_sf, p_symm, p_kv, wfc_basis, psi_in, stowf, nlpp_in);

for (int ipol = 0; ipol < 3; ++ipol)
{
Expand All @@ -65,7 +68,7 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,

if (ModuleSymmetry::Symmetry::symm_flag == 1)
{
p_symm->symmetrize_mat3(sigmatot, GlobalC::ucell.lat);
p_symm->symmetrize_mat3(sigmatot, this->ucell->lat);
}

bool ry = false;
Expand Down Expand Up @@ -108,7 +111,7 @@ void Sto_Stress_PW::sto_stress_kin(ModuleBase::matrix& sigma,
gk[0] = new double[npwx];
gk[1] = new double[npwx];
gk[2] = new double[npwx];
double tpiba = ModuleBase::TWO_PI / GlobalC::ucell.lat0;
double tpiba = ModuleBase::TWO_PI / this->ucell->lat0;
double twobysqrtpi = 2.0 / std::sqrt(ModuleBase::PI);
double* kfac = new double[npwx];
int nksbands = psi_in->get_nbands();
Expand Down Expand Up @@ -180,7 +183,7 @@ void Sto_Stress_PW::sto_stress_kin(ModuleBase::matrix& sigma,
{
for (int m = 0; m < 3; ++m)
{
s_kin(l, m) *= ModuleBase::e2 / GlobalC::ucell.omega;
s_kin(l, m) *= ModuleBase::e2 / this->ucell->omega;
}
}

Expand All @@ -196,7 +199,7 @@ void Sto_Stress_PW::sto_stress_kin(ModuleBase::matrix& sigma,
// do symmetry
if (ModuleSymmetry::Symmetry::symm_flag == 1)
{
p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
p_symm->symmetrize_mat3(sigma, this->ucell->lat);
}
delete[] gk[0];
delete[] gk[1];
Expand All @@ -215,12 +218,14 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
K_Vectors* p_kv,
ModulePW::PW_Basis_K* wfc_basis,
const psi::Psi<complex<double>>* psi_in,
Stochastic_WF& stowf)
Stochastic_WF& stowf,
pseudopot_cell_vnl* nlpp_in)
{
ModuleBase::TITLE("Sto_Stress_Func", "stres_nl");
ModuleBase::timer::tick("Sto_Stress_Func", "stres_nl");

const int nkb = GlobalC::ppcell.nkb;
this->nlpp = nlpp_in;
const int nkb = this->nlpp->nkb;
if (nkb == 0)
{
ModuleBase::timer::tick("Stress_Func", "stres_nl");
Expand Down Expand Up @@ -255,9 +260,9 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
if (GlobalV::NSPIN == 2)
GlobalV::CURRENT_SPIN = p_kv->isk[ik];
// generate vkb
if (GlobalC::ppcell.nkb > 0)
if (this->nlpp->nkb > 0)
{
GlobalC::ppcell.getvnl(ik, GlobalC::ppcell.vkb);
this->nlpp->getvnl(ik, this->nlpp->vkb);
}

// get becp according to wave functions and vkb
Expand All @@ -277,7 +282,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
&npmks,
&npw,
&ModuleBase::ONE,
GlobalC::ppcell.vkb.c,
this->nlpp->vkb.c,
&npwx,
psi_in->get_pointer(),
&npwx,
Expand All @@ -292,7 +297,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
&npmsto,
&npw,
&ModuleBase::ONE,
GlobalC::ppcell.vkb.c,
this->nlpp->vkb.c,
&npwx,
stowf.shchi->get_pointer(),
&npwx,
Expand Down Expand Up @@ -347,7 +352,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
// second termi
if (ipol == jpol)
{
pvkb = &GlobalC::ppcell.vkb(i, 0);
pvkb = &this->nlpp->vkb(i, 0);
for (int ig = 0; ig < npw; ++ig)
{
pdbecp_noevc[ig] -= pvkb[ig];
Expand All @@ -364,7 +369,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
else
qm1 = 0;
pdbecp_noevc[ig]
-= 2.0 * pvkb[ig] * qvec0[ipol][0] * qvec0[jpol][0] * qm1 * GlobalC::ucell.tpiba;
-= 2.0 * pvkb[ig] * qvec0[ipol][0] * qvec0[jpol][0] * qm1 * this->ucell->tpiba;
} // end ig
} // end i

Expand Down Expand Up @@ -406,14 +411,14 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
fac = p_kv->wk[ik];
int iat = 0;
int sum = 0;
for (int it = 0; it < GlobalC::ucell.ntype; ++it)
for (int it = 0; it < this->ucell->ntype; ++it)
{
const int Nprojs = GlobalC::ucell.atoms[it].ncpp.nh;
for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ++ia)
const int Nprojs = this->ucell->atoms[it].ncpp.nh;
for (int ia = 0; ia < this->ucell->atoms[it].na; ++ia)
{
for (int ip = 0; ip < Nprojs; ++ip)
{
double ps = GlobalC::ppcell.deeq(GlobalV::CURRENT_SPIN, iat, ip, ip);
double ps = this->nlpp->deeq(GlobalV::CURRENT_SPIN, iat, ip, ip);
const int inkb = sum + ip;
// out<<"\n ps = "<<ps;

Expand Down Expand Up @@ -447,7 +452,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
{
for (int jpol = 0; jpol < 3; ++jpol)
{
sigmanlc(ipol, jpol) *= 1.0 / GlobalC::ucell.omega;
sigmanlc(ipol, jpol) *= 1.0 / this->ucell->omega;
}
}

Expand All @@ -461,7 +466,7 @@ void Sto_Stress_PW::sto_stress_nl(ModuleBase::matrix& sigma,
// do symmetry
if (ModuleSymmetry::Symmetry::symm_flag == 1)
{
p_symm->symmetrize_mat3(sigma, GlobalC::ucell.lat);
p_symm->symmetrize_mat3(sigma, this->ucell->lat);
}

// this->print(ofs_running, "nonlocal stress", stresnl);
Expand Down
7 changes: 5 additions & 2 deletions source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,9 @@ class Sto_Stress_PW : public Stress_Func<double>
ModulePW::PW_Basis_K* wfc_basis,
const psi::Psi<complex<double>>* psi_in,
Stochastic_WF& stowf,
const Charge* const chr);
const Charge* const chr,
pseudopot_cell_vnl* nlpp_in,
const UnitCell& ucell_in);

private:
void sto_stress_kin(ModuleBase::matrix& sigma,
Expand All @@ -41,6 +43,7 @@ class Sto_Stress_PW : public Stress_Func<double>
K_Vectors* p_kv,
ModulePW::PW_Basis_K* wfc_basis,
const psi::Psi<complex<double>>* psi_in,
Stochastic_WF& stowf);
Stochastic_WF& stowf,
pseudopot_cell_vnl* nlpp_in);
};
#endif

0 comments on commit 9e9d748

Please sign in to comment.