Skip to content

Commit 6063015

Browse files
dyzhengdyzheng
committed
Fix: berry phase calculation in LCAO base with nscf calculation
1 parent 5f510c9 commit 6063015

File tree

8 files changed

+55
-35
lines changed

8 files changed

+55
-35
lines changed

‎source/input_conv.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -443,6 +443,10 @@ void Input_Conv::Convert(void)
443443
hsolver::HSolverLCAO::out_mat_hs = INPUT.out_mat_hs;
444444
hsolver::HSolverLCAO::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16
445445
elecstate::ElecStateLCAO::out_wfc_lcao = INPUT.out_wfc_lcao;
446+
if(INPUT.calculation == "nscf" && !INPUT.towannier90 && !INPUT.berry_phase)
447+
{
448+
elecstate::ElecStateLCAO::need_psi_grid = false;
449+
}
446450
#endif
447451

448452
GlobalC::en.dos_emin_ev = INPUT.dos_emin_ev;

‎source/module_elecstate/elecstate_lcao.cpp‎

Lines changed: 31 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
namespace elecstate
99
{
1010
int ElecStateLCAO::out_wfc_lcao = 0;
11+
bool ElecStateLCAO::need_psi_grid = 1;
1112

1213
// multi-k case
1314
void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
@@ -38,25 +39,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
3839
for (int ik = 0; ik < psi.get_nk(); ik++)
3940
{
4041
psi.fix_k(ik);
41-
this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), this->lowf->wfc_k_grid[ik], ik, this->ekb, this->wg);
42-
//added by zhengdy-soc, rearrange the wfc_k_grid from [up,down,up,down...] to [up,up...down,down...],
43-
if(GlobalV::NSPIN==4)
44-
{
45-
int row = GlobalC::GridT.lgd;
46-
std::vector<std::complex<double>> tmp(row);
47-
for(int ib=0; ib<GlobalV::NBANDS; ib++)
48-
{
49-
for(int iw=0; iw<row / GlobalV::NPOL; iw++)
50-
{
51-
tmp[iw] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL];
52-
tmp[iw + row / GlobalV::NPOL] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL + 1];
53-
}
54-
for(int iw=0; iw<row; iw++)
55-
{
56-
this->lowf->wfc_k_grid[ik][ib][iw] = tmp[iw];
57-
}
58-
}
59-
}
42+
this->print_psi(psi);
6043
}
6144
}
6245

@@ -112,8 +95,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
11295
if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx")
11396
{
11497
psi.fix_k(ik);
115-
double** wfc_grid = nullptr; // output but not do "2d-to-grid" conversion
116-
this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), wfc_grid, this->ekb, this->wg);
98+
this->print_psi(psi);
11799
}
118100
//this->loc->dm2dToGrid(this->loc->dm_gamma[ik], this->loc->DM[ik]); // transform dm_gamma[is].c to this->loc->DM[is]
119101
this->loc->cal_dk_gamma_from_2D_pub();
@@ -158,11 +140,36 @@ void ElecStateLCAO::print_psi(const psi::Psi<double>& psi_in)
158140
}
159141
void ElecStateLCAO::print_psi(const psi::Psi<std::complex<double>>& psi_in)
160142
{
161-
if(!ElecStateLCAO::out_wfc_lcao) return;
143+
if(!ElecStateLCAO::out_wfc_lcao && !ElecStateLCAO::need_psi_grid) return;
162144

163145
// output but not do "2d-to-grid" conversion
164-
std::complex<double>** wfc_grid = nullptr;
165-
this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi_in.get_pointer(), wfc_grid, psi_in.get_current_k(), this->ekb, this->wg);
146+
std::complex<double>** wfc_grid = nullptr;
147+
int ik = psi_in.get_current_k();
148+
if(ElecStateLCAO::need_psi_grid)
149+
{
150+
wfc_grid = this->lowf->wfc_k_grid[ik];
151+
}
152+
this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi_in.get_pointer(), wfc_grid, ik, this->ekb, this->wg);
153+
154+
//added by zhengdy-soc, rearrange the wfc_k_grid from [up,down,up,down...] to [up,up...down,down...],
155+
if(ElecStateLCAO::need_psi_grid && GlobalV::NSPIN==4)
156+
{
157+
int row = GlobalC::GridT.lgd;
158+
std::vector<std::complex<double>> tmp(row);
159+
for(int ib=0; ib<GlobalV::NBANDS; ib++)
160+
{
161+
for(int iw=0; iw<row / GlobalV::NPOL; iw++)
162+
{
163+
tmp[iw] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL];
164+
tmp[iw + row / GlobalV::NPOL] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL + 1];
165+
}
166+
for(int iw=0; iw<row; iw++)
167+
{
168+
this->lowf->wfc_k_grid[ik][ib][iw] = tmp[iw];
169+
}
170+
}
171+
}
172+
166173
return;
167174
}
168175

‎source/module_elecstate/elecstate_lcao.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@ class ElecStateLCAO : public ElecState
4141
virtual void print_psi(const psi::Psi<std::complex<double>>& psi_in)override;
4242

4343
static int out_wfc_lcao;
44+
static bool need_psi_grid;
4445

4546
private:
4647
// calculate electronic charge density on grid points or density matrix in real space

‎source/module_esolver/esolver_ks_lcao_elec.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -395,7 +395,7 @@ namespace ModuleESolver
395395
if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag == 0)
396396
{
397397
berryphase bp(this->LOWF);
398-
bp.Macroscopic_polarization(nullptr);
398+
bp.Macroscopic_polarization(this->psi);
399399
}
400400

401401
return;

‎source/src_io/berryphase.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -325,7 +325,7 @@ double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::co
325325
if(GlobalV::NSPIN!=4)
326326
{
327327
//std::complex<double> my_det = lcao_method.det_berryphase(ik_1,ik_2,dk,nbands);
328-
zeta = zeta * lcao_method.det_berryphase(ik_1,ik_2,dk,nbands, *this->lowf);
328+
zeta = zeta * lcao_method.det_berryphase(ik_1,ik_2,dk,nbands, *this->lowf, psi_in);
329329
// test by jingan
330330
//GlobalV::ofs_running << "methon 1: det = " << my_det << std::endl;
331331
// test by jingan
@@ -456,7 +456,7 @@ void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>*
456456
GlobalV::ofs_running << "\n\n\n\n";
457457
GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
458458
GlobalV::ofs_running << " | |" << std::endl;
459-
GlobalV::ofs_running << " | POLARIZATION GlobalV::CALCULATION: |" << std::endl;
459+
GlobalV::ofs_running << " | POLARIZATION CALCULATION: |" << std::endl;
460460
GlobalV::ofs_running << " | Modern Theory of Polarization |" << std::endl;
461461
GlobalV::ofs_running << " | calculate the Macroscopic polarization of a crystalline insulator |" << std::endl;
462462
GlobalV::ofs_running << " | by using Berry Phase method. |" << std::endl;

‎source/src_io/unk_overlap_lcao.cpp‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -723,7 +723,8 @@ void unkOverlap_lcao::prepare_midmatrix_pblas(const int ik_L, const int ik_R, co
723723

724724
std::complex<double> unkOverlap_lcao::det_berryphase(const int ik_L, const int ik_R,
725725
const ModuleBase::Vector3<double> dk, const int occ_bands,
726-
Local_Orbital_wfc &lowf)
726+
Local_Orbital_wfc &lowf,
727+
const psi::Psi<std::complex<double>>* psi_in)
727728
{
728729
const std::complex<double> minus = std::complex<double>(-1.0,0.0);
729730
std::complex<double> det = std::complex<double>(1.0,0.0);
@@ -743,14 +744,14 @@ std::complex<double> unkOverlap_lcao::det_berryphase(const int ik_L, const int i
743744
int one = 1;
744745
#ifdef __MPI
745746
pzgemm_(&transa,&transb,&occBands,&nlocal,&nlocal,&alpha[0],
746-
lowf.wfc_k.at(ik_L).c,&one,&one,lowf.ParaV->desc,
747+
&psi_in[0](ik_L, 0, 0), &one, &one, lowf.ParaV->desc,
747748
midmatrix,&one,&one,lowf.ParaV->desc,
748749
&beta[0],
749750
C_matrix,&one,&one,lowf.ParaV->desc);
750751

751752
pzgemm_(&transb,&transb,&occBands,&occBands,&nlocal,&alpha[0],
752753
C_matrix,&one,&one,lowf.ParaV->desc,
753-
lowf.wfc_k.at(ik_R).c,&one,&one,lowf.ParaV->desc,
754+
&psi_in[0](ik_R, 0, 0), &one, &one, lowf.ParaV->desc,
754755
&beta[0],
755756
out_matrix,&one,&one,lowf.ParaV->desc);
756757

‎source/src_io/unk_overlap_lcao.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,8 @@ class unkOverlap_lcao
6868
void prepare_midmatrix_pblas(const int ik_L, const int ik_R, const ModuleBase::Vector3<double> dk, std::complex<double> *&midmatrix, const Parallel_Orbitals &pv);
6969
std::complex<double> det_berryphase(const int ik_L, const int ik_R,
7070
const ModuleBase::Vector3<double> dk, const int occ_bands,
71-
Local_Orbital_wfc &lowf);
71+
Local_Orbital_wfc &lowf,
72+
const psi::Psi<std::complex<double>>* psi_in);
7273

7374
void test(std::complex<double>*** wfc_k_grid);
7475

‎source/src_pw/klist.cpp‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include "../src_parallel/parallel_reduce.h"
66
#include "../src_parallel/parallel_common.h"
77
#include "../module_base/memory.h"
8+
#include "src_io/berryphase.h"
89

910
K_Vectors::K_Vectors()
1011
{
@@ -80,11 +81,16 @@ void K_Vectors::set(
8081
}
8182

8283
// (2)
83-
this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag);
84-
if(ModuleSymmetry::Symmetry::symm_flag || is_mp)
84+
//only berry phase need all kpoints including time-reversal symmetry!
85+
//if symm_flag is not set, only time-reversal symmetry would be considered.
86+
if(!berryphase::berry_phase_flag)
8587
{
86-
this->update_use_ibz();
87-
this->nks = this->nkstot = this->nkstot_ibz;
88+
this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag);
89+
if(ModuleSymmetry::Symmetry::symm_flag || is_mp)
90+
{
91+
this->update_use_ibz();
92+
this->nks = this->nkstot = this->nkstot_ibz;
93+
}
8894
}
8995

9096
// (3)

0 commit comments

Comments
 (0)