Skip to content

Commit 3e22ff9

Browse files
dyzhengdyzheng
committed
Refactor: add elecstate.cpp and elecstate_pw.cpp, HSolverPW is finished
1 parent 7e42ed7 commit 3e22ff9

File tree

15 files changed

+431
-74
lines changed

15 files changed

+431
-74
lines changed

‎source/module_elecstate/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,5 @@ add_library(
22
elecstate
33
OBJECT
44
elecstate_pw.cpp
5+
elecstate.cpp
56
)

‎source/module_elecstate/elecstate.cpp‎

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
#include "elecstate.h"
2+
#include "src_pw/occupy.h"
3+
#include "module_base/global_variable.h"
4+
#include "module_base/tool_title.h"
5+
#include "src_pw/global.h"
6+
#include "src_parallel/parallel_reduce.h"
7+
8+
namespace elecstate
9+
{
10+
11+
const double* ElecState::getRho(int spin) const
12+
{
13+
//hamilt::MatrixBlock<double> temp{&(this->charge->rho[spin][0]), 1, this->charge->nrxx}; // this->chr->get_nspin(), this->chr->get_nrxx()};
14+
return &(this->charge->rho[spin][0]);
15+
}
16+
17+
void ElecState::calculate_weights(void)
18+
{
19+
ModuleBase::TITLE("ElecState", "calculate_weights");
20+
21+
// for test
22+
// std::cout << " gaussian_broadening = " << use_gaussian_broadening << std::endl;
23+
// std::cout << " tetrahedron_method = " << use_tetrahedron_method << std::endl;
24+
// std::cout << " fixed_occupations = " << fixed_occupations << std::endl;
25+
double** ekb_tmp = new double*[this->ekb.nr];
26+
for(int i=0; i<this->ekb.nr; ++i)
27+
{
28+
ekb_tmp[i] = &(this->ekb(i, 0));
29+
}
30+
31+
if (GlobalV::KS_SOLVER == "selinv")
32+
{
33+
GlobalV::ofs_running << " Could not calculate occupation." << std::endl;
34+
return;
35+
}
36+
37+
if (!Occupy::use_gaussian_broadening && !Occupy::use_tetrahedron_method && !Occupy::fixed_occupations)
38+
{
39+
if (GlobalV::TWO_EFERMI)
40+
{
41+
Occupy::iweights(GlobalC::kv.nks,
42+
GlobalC::kv.wk,
43+
GlobalV::NBANDS,
44+
GlobalC::ucell.magnet.get_nelup(),
45+
ekb_tmp,
46+
GlobalC::en.ef_up,
47+
this->wg,
48+
0,
49+
GlobalC::kv.isk);
50+
Occupy::iweights(GlobalC::kv.nks,
51+
GlobalC::kv.wk,
52+
GlobalV::NBANDS,
53+
GlobalC::ucell.magnet.get_neldw(),
54+
ekb_tmp,
55+
GlobalC::en.ef_dw,
56+
this->wg,
57+
1,
58+
GlobalC::kv.isk);
59+
// ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
60+
}
61+
else
62+
{
63+
// -1 means don't need to consider spin.
64+
Occupy::iweights(GlobalC::kv.nks,
65+
GlobalC::kv.wk,
66+
GlobalV::NBANDS,
67+
GlobalC::CHR.nelec,
68+
ekb_tmp,
69+
this->ef,
70+
this->wg,
71+
-1,
72+
GlobalC::kv.isk);
73+
}
74+
}
75+
else if (Occupy::use_tetrahedron_method)
76+
{
77+
ModuleBase::WARNING_QUIT("calculate_weights", "not implemented yet,coming soon!");
78+
// if(my_rank == 0)
79+
// {
80+
// tweights(GlobalC::kv.nkstot, nspin, GlobalV::NBANDS, GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
81+
//this->ef, this->wg);
82+
// }
83+
}
84+
else if (Occupy::use_gaussian_broadening)
85+
{
86+
if (GlobalV::TWO_EFERMI)
87+
{
88+
double demet_up = 0.0;
89+
double demet_dw = 0.0;
90+
Occupy::gweights(GlobalC::kv.nks,
91+
GlobalC::kv.wk,
92+
GlobalV::NBANDS,
93+
GlobalC::ucell.magnet.get_nelup(),
94+
Occupy::gaussian_parameter,
95+
Occupy::gaussian_type,
96+
ekb_tmp,
97+
GlobalC::en.ef_up,
98+
demet_up,
99+
this->wg,
100+
0,
101+
GlobalC::kv.isk);
102+
Occupy::gweights(GlobalC::kv.nks,
103+
GlobalC::kv.wk,
104+
GlobalV::NBANDS,
105+
GlobalC::ucell.magnet.get_neldw(),
106+
Occupy::gaussian_parameter,
107+
Occupy::gaussian_type,
108+
ekb_tmp,
109+
GlobalC::en.ef_dw,
110+
demet_dw,
111+
this->wg,
112+
1,
113+
GlobalC::kv.isk);
114+
GlobalC::en.demet = demet_up + demet_dw;
115+
}
116+
else
117+
{
118+
// -1 means is no related to spin.
119+
Occupy::gweights(GlobalC::kv.nks,
120+
GlobalC::kv.wk,
121+
GlobalV::NBANDS,
122+
GlobalC::CHR.nelec,
123+
Occupy::gaussian_parameter,
124+
Occupy::gaussian_type,
125+
ekb_tmp,
126+
this->ef,
127+
GlobalC::en.demet,
128+
this->wg,
129+
-1,
130+
GlobalC::kv.isk);
131+
}
132+
133+
// qianrui fix a bug on 2021-7-21
134+
Parallel_Reduce::reduce_double_allpool(GlobalC::en.demet);
135+
}
136+
else if (Occupy::fixed_occupations)
137+
{
138+
// fix occupations need nelup and neldw.
139+
// mohan add 2011-04-03
140+
this->ef = -1.0e+20;
141+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
142+
{
143+
for (int ibnd = 0; ibnd < GlobalV::NBANDS; ibnd++)
144+
{
145+
if (this->wg(ik, ibnd) > 0.0)
146+
{
147+
this->ef = std::max(this->ef, ekb_tmp[ik][ibnd]);
148+
}
149+
}
150+
}
151+
}
152+
153+
if (GlobalV::TWO_EFERMI)
154+
{
155+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_up);
156+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_dw);
157+
}
158+
else
159+
{
160+
double ebotom = ekb_tmp[0][0];
161+
double etop = ekb_tmp[0][0];
162+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
163+
{
164+
for (int ib = 0; ib < GlobalV::NBANDS; ib++)
165+
{
166+
ebotom = min(ebotom, ekb_tmp[ik][ib]);
167+
etop = max(etop, ekb_tmp[ik][ib]);
168+
}
169+
}
170+
171+
// parallel
172+
Parallel_Reduce::gather_max_double_all(this->ef);
173+
Parallel_Reduce::gather_max_double_all(etop);
174+
Parallel_Reduce::gather_min_double_all(ebotom);
175+
176+
// not parallel yet!
177+
// OUT(GlobalV::ofs_running,"Top Energy (eV)", etop * ModuleBase::Ry_to_eV);
178+
// OUT(GlobalV::ofs_running,"Fermi Energy (eV)", this->ef * ModuleBase::Ry_to_eV);
179+
// OUT(GlobalV::ofs_running,"Bottom Energy (eV)", ebotom * ModuleBase::Ry_to_eV);
180+
// OUT(GlobalV::ofs_running,"Range Energy (eV)", etop-ebotom * ModuleBase::Ry_to_eV);
181+
}
182+
183+
delete[] ekb_tmp;
184+
185+
return;
186+
}
187+
188+
} // namespace elecstate

‎source/module_elecstate/elecstate.h‎

Lines changed: 29 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,28 +11,47 @@ namespace elecstate
1111
class ElecState
1212
{
1313
public:
14-
virtual void init(Charge *chg_in
15-
/*const Basis &basis, const Cell &cell*/)
16-
= 0;
14+
virtual void init(
15+
Charge *chg_in, //pointer for class Charge
16+
int nk_in, //number of k points
17+
int nb_in) //number of bands
18+
{
19+
this->charge = chg_in;
20+
this->ekb.create(nk_in, nb_in);
21+
this->wg.create(nk_in, nb_in);
22+
}
1723

1824
// return current electronic density rho, as a input for constructing Hamiltonian
19-
virtual const hamilt::MatrixBlock<double> getRho() const = 0;
25+
virtual const double* getRho(int spin) const;
2026

2127
// calculate electronic charge density on grid points or density matrix in real space
2228
// the consequence charge density rho saved into rho_out, preparing for charge mixing.
23-
virtual void updateRhoK(const psi::Psi<std::complex<double>> &psi) = 0;
24-
virtual void updateRhoK(const psi::Psi<double> &psi)
25-
{
26-
return;
27-
}
29+
virtual void psiToRho(const psi::Psi<std::complex<double>> &psi){ return; }
30+
virtual void psiToRho(const psi::Psi<double> &psi){ return; }
31+
//virtual void updateRhoK(const psi::Psi<std::complex<double>> &psi) = 0;
32+
//virtual void updateRhoK(const psi::Psi<double> &psi)=0
2833

2934
// update charge density for next scf step
3035
// in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge
3136
// density rho among these rho,
3237
// 1. input rho would be store to file for restart
3338
// 2. calculated rho should be near with input rho when convergence has achieved
3439
// 3. new rho should be input rho for next scf step.
35-
virtual void getNewRho() = 0;
40+
virtual void getNewRho(){return;}
41+
42+
//calculate wg from ekb
43+
void calculate_weights(void);
44+
45+
Charge* charge;
46+
//energy for sum of electrons
47+
double eband = 0.0;
48+
//Fermi energy
49+
double ef = 0.0;
50+
51+
// band energy at each k point, each band.
52+
ModuleBase::matrix ekb;
53+
// occupation weight for each k-point and band
54+
ModuleBase::matrix wg;
3655
};
3756

3857
} // namespace elecstate

0 commit comments

Comments
 (0)