Skip to content

Commit 8f15d2a

Browse files
pxlxingliangpxlxingliang
committed
test(hsolver): add the test of ELPA/SCALAPACK
1 parent af32cfe commit 8f15d2a

27 files changed

+3484
-80
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: 250 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,250 @@
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+
int nbands = this->ekb.nc;
31+
32+
if (GlobalV::KS_SOLVER == "selinv")
33+
{
34+
GlobalV::ofs_running << " Could not calculate occupation." << std::endl;
35+
return;
36+
}
37+
38+
if (!Occupy::use_gaussian_broadening && !Occupy::use_tetrahedron_method && !Occupy::fixed_occupations)
39+
{
40+
if (GlobalV::TWO_EFERMI)
41+
{
42+
Occupy::iweights(GlobalC::kv.nks,
43+
GlobalC::kv.wk,
44+
nbands,
45+
GlobalC::ucell.magnet.get_nelup(),
46+
ekb_tmp,
47+
GlobalC::en.ef_up,
48+
this->wg,
49+
0,
50+
GlobalC::kv.isk);
51+
Occupy::iweights(GlobalC::kv.nks,
52+
GlobalC::kv.wk,
53+
nbands,
54+
GlobalC::ucell.magnet.get_neldw(),
55+
ekb_tmp,
56+
GlobalC::en.ef_dw,
57+
this->wg,
58+
1,
59+
GlobalC::kv.isk);
60+
// ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
61+
}
62+
else
63+
{
64+
// -1 means don't need to consider spin.
65+
Occupy::iweights(GlobalC::kv.nks,
66+
GlobalC::kv.wk,
67+
nbands,
68+
GlobalC::CHR.nelec,
69+
ekb_tmp,
70+
this->ef,
71+
this->wg,
72+
-1,
73+
GlobalC::kv.isk);
74+
}
75+
}
76+
else if (Occupy::use_tetrahedron_method)
77+
{
78+
ModuleBase::WARNING_QUIT("calculate_weights", "not implemented yet,coming soon!");
79+
// if(my_rank == 0)
80+
// {
81+
// tweights(GlobalC::kv.nkstot, nspin, nbands, GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
82+
//this->ef, this->wg);
83+
// }
84+
}
85+
else if (Occupy::use_gaussian_broadening)
86+
{
87+
if (GlobalV::TWO_EFERMI)
88+
{
89+
double demet_up = 0.0;
90+
double demet_dw = 0.0;
91+
Occupy::gweights(GlobalC::kv.nks,
92+
GlobalC::kv.wk,
93+
nbands,
94+
GlobalC::ucell.magnet.get_nelup(),
95+
Occupy::gaussian_parameter,
96+
Occupy::gaussian_type,
97+
ekb_tmp,
98+
GlobalC::en.ef_up,
99+
demet_up,
100+
this->wg,
101+
0,
102+
GlobalC::kv.isk);
103+
Occupy::gweights(GlobalC::kv.nks,
104+
GlobalC::kv.wk,
105+
nbands,
106+
GlobalC::ucell.magnet.get_neldw(),
107+
Occupy::gaussian_parameter,
108+
Occupy::gaussian_type,
109+
ekb_tmp,
110+
GlobalC::en.ef_dw,
111+
demet_dw,
112+
this->wg,
113+
1,
114+
GlobalC::kv.isk);
115+
GlobalC::en.demet = demet_up + demet_dw;
116+
}
117+
else
118+
{
119+
// -1 means is no related to spin.
120+
Occupy::gweights(GlobalC::kv.nks,
121+
GlobalC::kv.wk,
122+
nbands,
123+
GlobalC::CHR.nelec,
124+
Occupy::gaussian_parameter,
125+
Occupy::gaussian_type,
126+
ekb_tmp,
127+
this->ef,
128+
GlobalC::en.demet,
129+
this->wg,
130+
-1,
131+
GlobalC::kv.isk);
132+
}
133+
134+
// qianrui fix a bug on 2021-7-21
135+
Parallel_Reduce::reduce_double_allpool(GlobalC::en.demet);
136+
}
137+
else if (Occupy::fixed_occupations)
138+
{
139+
// fix occupations need nelup and neldw.
140+
// mohan add 2011-04-03
141+
this->ef = -1.0e+20;
142+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
143+
{
144+
for (int ibnd = 0; ibnd < nbands; ibnd++)
145+
{
146+
if (this->wg(ik, ibnd) > 0.0)
147+
{
148+
this->ef = std::max(this->ef, ekb_tmp[ik][ibnd]);
149+
}
150+
}
151+
}
152+
}
153+
154+
if (GlobalV::TWO_EFERMI)
155+
{
156+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_up);
157+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_dw);
158+
}
159+
else
160+
{
161+
double ebotom = ekb_tmp[0][0];
162+
double etop = ekb_tmp[0][0];
163+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
164+
{
165+
for (int ib = 0; ib < nbands; ib++)
166+
{
167+
ebotom = min(ebotom, ekb_tmp[ik][ib]);
168+
etop = max(etop, ekb_tmp[ik][ib]);
169+
}
170+
}
171+
172+
// parallel
173+
Parallel_Reduce::gather_max_double_all(this->ef);
174+
Parallel_Reduce::gather_max_double_all(etop);
175+
Parallel_Reduce::gather_min_double_all(ebotom);
176+
177+
// not parallel yet!
178+
// OUT(GlobalV::ofs_running,"Top Energy (eV)", etop * ModuleBase::Ry_to_eV);
179+
// OUT(GlobalV::ofs_running,"Fermi Energy (eV)", this->ef * ModuleBase::Ry_to_eV);
180+
// OUT(GlobalV::ofs_running,"Bottom Energy (eV)", ebotom * ModuleBase::Ry_to_eV);
181+
// OUT(GlobalV::ofs_running,"Range Energy (eV)", etop-ebotom * ModuleBase::Ry_to_eV);
182+
}
183+
184+
delete[] ekb_tmp;
185+
186+
return;
187+
}
188+
189+
double ElecState::eBandK(const int& ik)
190+
{
191+
ModuleBase::TITLE("ElecStatePW","eBandK");
192+
double eband_k = 0.0;
193+
for (int ibnd = 0;ibnd < this->ekb.nc; ibnd++)
194+
{
195+
eband_k += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
196+
}
197+
return eband_k;
198+
}
199+
200+
void ElecState::print_band(const int &ik, const int &printe, const int &iter)
201+
{
202+
//check the band energy.
203+
bool wrong = false;
204+
int nbands = this->ekb.nc;
205+
for(int ib=0; ib<nbands; ++ib)
206+
{
207+
if( abs( this->ekb(ik, ib) ) > 1.0e10)
208+
{
209+
GlobalV::ofs_warning << " ik=" << ik+1 << " ib=" << ib+1 << " " << this->ekb(ik, ib) << " Ry" << std::endl;
210+
wrong = true;
211+
}
212+
}
213+
if(wrong)
214+
{
215+
ModuleBase::WARNING_QUIT("Threshold_Elec::print_eigenvalue","Eigenvalues are too large!");
216+
}
217+
218+
219+
220+
if(GlobalV::MY_RANK==0)
221+
{
222+
//if( GlobalV::DIAGO_TYPE == "selinv" ) xiaohui modify 2013-09-02
223+
if(GlobalV::KS_SOLVER=="selinv") //xiaohui add 2013-09-02
224+
{
225+
GlobalV::ofs_running << " No eigenvalues are available for selected inversion methods." << std::endl;
226+
}
227+
else
228+
{
229+
if( printe>0 && ((iter+1) % printe == 0))
230+
{
231+
// NEW_PART("ENERGY BANDS (Rydberg), (eV)");
232+
GlobalV::ofs_running << std::setprecision(6);
233+
GlobalV::ofs_running << " Energy (eV) & Occupations for spin=" << GlobalV::CURRENT_SPIN+1 << " K-point=" << ik+1 << std::endl;
234+
GlobalV::ofs_running << std::setiosflags(ios::showpoint);
235+
for(int ib=0;ib<nbands;ib++)
236+
{
237+
GlobalV::ofs_running << " "<< std::setw(6) << ib+1
238+
<< std::setw(15) << this->ekb(ik, ib) * ModuleBase::Ry_to_eV;
239+
// for the first electron iteration, we don't have the energy
240+
// spectrum, so we can't get the occupations.
241+
GlobalV::ofs_running << std::setw(15) << this->wg(ik,ib);
242+
GlobalV::ofs_running << std::endl;
243+
}
244+
}
245+
}
246+
}
247+
return;
248+
}
249+
250+
} // namespace elecstate

‎source/module_elecstate/elecstate.h‎

Lines changed: 36 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
#ifndef ELECSTATE_H
22
#define ELECSTATE_H
33

4-
#include "module_hamilt/matrixblock.h"
54
#include "module_psi/psi.h"
65
#include "src_pw/charge.h"
76

@@ -11,28 +10,54 @@ namespace elecstate
1110
class ElecState
1211
{
1312
public:
14-
virtual void init(Charge *chg_in
15-
/*const Basis &basis, const Cell &cell*/)
16-
= 0;
13+
virtual void init(
14+
Charge *chg_in, //pointer for class Charge
15+
int nk_in, //number of k points
16+
int nb_in) //number of bands
17+
{
18+
this->charge = chg_in;
19+
this->ekb.create(nk_in, nb_in);
20+
this->wg.create(nk_in, nb_in);
21+
}
1722

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

2126
// calculate electronic charge density on grid points or density matrix in real space
2227
// 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-
}
28+
virtual void psiToRho(const psi::Psi<std::complex<double>> &psi){ return; }
29+
virtual void psiToRho(const psi::Psi<double> &psi){ return; }
30+
//virtual void updateRhoK(const psi::Psi<std::complex<double>> &psi) = 0;
31+
//virtual void updateRhoK(const psi::Psi<double> &psi)=0
2832

2933
// update charge density for next scf step
3034
// in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge
3135
// density rho among these rho,
3236
// 1. input rho would be store to file for restart
3337
// 2. calculated rho should be near with input rho when convergence has achieved
3438
// 3. new rho should be input rho for next scf step.
35-
virtual void getNewRho() = 0;
39+
virtual void getNewRho(){return;}
40+
41+
//calculate wg from ekb
42+
void calculate_weights(void);
43+
44+
Charge* charge;
45+
//energy for sum of electrons
46+
double eband = 0.0;
47+
//Fermi energy
48+
double ef = 0.0;
49+
50+
// band energy at each k point, each band.
51+
ModuleBase::matrix ekb;
52+
// occupation weight for each k-point and band
53+
ModuleBase::matrix wg;
54+
55+
protected:
56+
// calculate ebands for each k point
57+
double eBandK(const int& ik);
58+
59+
//print and check for band energy and occupations
60+
void print_band(const int &ik, const int &printe, const int &iter);
3661
};
3762

3863
} // namespace elecstate

0 commit comments

Comments
 (0)