Skip to content

Commit

Permalink
Fix: wrong results of wannier function
Browse files Browse the repository at this point in the history
  • Loading branch information
@Qianruipku
Qianruipku committed Aug 30, 2023
1 parent e9975c8 commit 3554fc5
Show file tree
Hide file tree
Showing 3 changed files with 14 additions and 21 deletions.
1 change: 0 additions & 1 deletion source/module_esolver/esolver_ks_lcao_elec.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -568,7 +568,6 @@ namespace ModuleESolver
toWannier90 myWannier(this->kv.nkstot, GlobalC::ucell.G, this->LOWF.wfc_k_grid);
myWannier.init_wannier_lcao(this->GridT,
this->pelec->ekb,
this->pw_rho,
this->pw_wfc,
this->pw_big,
this->sf,
Expand Down
26 changes: 12 additions & 14 deletions source/module_io/to_wannier90.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ void toWannier90::init_wannier_pw(const ModuleBase::matrix& ekb,

writeUNK(wfcpw, *psi, bigpw);
outEIG(ekb);
cal_Mmn(*psi, rhopw, wfcpw);
cal_Mmn(*psi, wfcpw);
cal_Amn(*psi, wfcpw);

/*
Expand Down Expand Up @@ -93,7 +93,6 @@ void toWannier90::init_wannier_pw(const ModuleBase::matrix& ekb,
#ifdef __LCAO
void toWannier90::init_wannier_lcao(const Grid_Technique& gt,
const ModuleBase::matrix& ekb,
const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw,
const ModulePW::PW_Basis_Big* bigpw,
const Structure_Factor& sf,
Expand All @@ -118,7 +117,7 @@ void toWannier90::init_wannier_lcao(const Grid_Technique& gt,

getUnkFromLcao(wfcpw, sf, kv, wfcpw->npwk_max);
cal_Amn(this->unk_inLcao[0], wfcpw);
cal_Mmn(this->unk_inLcao[0], rhopw, wfcpw);
cal_Mmn(this->unk_inLcao[0], wfcpw);
writeUNK(wfcpw, this->unk_inLcao[0], bigpw);
outEIG(ekb);
}
Expand Down Expand Up @@ -699,7 +698,6 @@ void toWannier90::cal_Amn(const psi::Psi<std::complex<double>>& psi_pw, const Mo
}

void toWannier90::cal_Mmn(const psi::Psi<std::complex<double>>& psi_pw,
const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw)
{
// test by jingan
Expand Down Expand Up @@ -775,7 +773,7 @@ void toWannier90::cal_Mmn(const psi::Psi<std::complex<double>>& psi_pw,
// std::complex<double> *unk_L_r = new std::complex<double>[wfcpw->nrxx];
// ToRealSpace(cal_ik,n,psi_pw,unk_L_r,phase_G);
// mmn = unkdotb(unk_L_r,cal_ikb,m,psi_pw);
mmn = unkdotkb(rhopw, wfcpw, cal_ik, cal_ikb, n, m, phase_G, psi_pw);
mmn = unkdotkb(wfcpw, cal_ik, cal_ikb, n, m, phase_G, psi_pw);
// delete[] unk_L_r;
}
else
Expand Down Expand Up @@ -1622,8 +1620,7 @@ std::complex<double> toWannier90::unkdotb(const std::complex<double> *psir,
return result;
}
*/
std::complex<double> toWannier90::unkdotkb(const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw,
std::complex<double> toWannier90::unkdotkb(const ModulePW::PW_Basis_K* wfcpw,
const int& ik,
const int& ikb,
const int& iband_L,
Expand All @@ -1633,33 +1630,34 @@ std::complex<double> toWannier90::unkdotkb(const ModulePW::PW_Basis* rhopw,
{
// (1) set value
std::complex<double> result(0.0, 0.0);
std::complex<double> *psir = new std::complex<double>[wfcpw->nmaxgr];
std::complex<double>* phase = new std::complex<double>[rhopw->nmaxgr];
std::complex<double>* psir = new std::complex<double>[wfcpw->nmaxgr];
std::complex<double>* phase = new std::complex<double>[wfcpw->nmaxgr];
ModuleBase::GlobalFunc::ZEROS(phase, wfcpw->nmaxgr);

// get the phase value in realspace
for (int ig = 0; ig < rhopw->npw; ig++)
for (int ig = 0; ig < wfcpw->npwk[ik]; ig++)
{
if (rhopw->gdirect[ig] == G) // It should be used carefully. We cannot judge if two double are equal.
if (wfcpw->getgdirect(ik,ig) == G) // It should be used carefully. We cannot judge if two double are equal.
{
phase[ig] = std::complex<double>(1.0, 0.0);
break;
}
}

// (2) fft and get value
rhopw->recip2real(phase, phase);
wfcpw->recip2real(phase, phase, ik);
wfcpw->recip2real(&psi_pw(ik, iband_L, 0), psir, ik);

for (int ir = 0; ir < wfcpw->nrxx; ir++)
{
psir[ir] *= phase[ir];
}

wfcpw->real2recip(psir, psir, ik);
wfcpw->real2recip(psir, psir, ikb); //ikb, not ik

std::complex<double> result_tem(0.0, 0.0);

for (int ig = 0; ig < psi_pw.get_ngk(ikb); ig++)
for (int ig = 0; ig <wfcpw->npwk[ikb]; ig++)
{
result_tem = result_tem + conj(psir[ig]) * psi_pw(ikb, iband_R, ig);
}
Expand Down
8 changes: 2 additions & 6 deletions source/module_io/to_wannier90.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,6 @@ class toWannier90
#ifdef __LCAO
void init_wannier_lcao(const Grid_Technique& gt,
const ModuleBase::matrix& ekb,
const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw,
const ModulePW::PW_Basis_Big* bigpw,
const Structure_Factor& sf,
Expand All @@ -89,9 +88,7 @@ class toWannier90
void read_nnkp(const K_Vectors& kv);
void outEIG(const ModuleBase::matrix& ekb);
void cal_Amn(const psi::Psi<std::complex<double>>& psi_pw, const ModulePW::PW_Basis_K* wfcpw);
void cal_Mmn(const psi::Psi<std::complex<double>>& psi_pw,
const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw);
void cal_Mmn(const psi::Psi<std::complex<double>>& psi_pw, const ModulePW::PW_Basis_K* wfcpw);
void produce_trial_in_pw(const psi::Psi<std::complex<double>>& psi_pw,
const int& ik,
const ModulePW::PW_Basis_K* wfcpw,
Expand All @@ -115,8 +112,7 @@ class toWannier90
// void ToRealSpace(const int &ik, const int &ib, const ModuleBase::ComplexMatrix *evc, std::complex<double> *psir,
// const ModuleBase::Vector3<double> G); std::complex<double> unkdotb(const std::complex<double> *psir, const int
// ikb, const int bandindex, const ModuleBase::ComplexMatrix *psi_pw);
std::complex<double> unkdotkb(const ModulePW::PW_Basis* rhopw,
const ModulePW::PW_Basis_K* wfcpw,
std::complex<double> unkdotkb(const ModulePW::PW_Basis_K* wfcpw,
const int& ik,
const int& ikb,
const int& iband_L,
Expand Down

0 comments on commit 3554fc5

Please sign in to comment.