@@ -27,6 +27,7 @@ void ElecState::calculate_weights(void)
2727 {
2828 ekb_tmp[i] = &(this ->ekb (i, 0 ));
2929 }
30+ int nbands = this ->ekb .nc ;
3031
3132 if (GlobalV::KS_SOLVER == " selinv" )
3233 {
@@ -40,7 +41,7 @@ void ElecState::calculate_weights(void)
4041 {
4142 Occupy::iweights (GlobalC::kv.nks ,
4243 GlobalC::kv.wk ,
43- GlobalV::NBANDS ,
44+ nbands ,
4445 GlobalC::ucell.magnet .get_nelup (),
4546 ekb_tmp,
4647 GlobalC::en.ef_up ,
@@ -49,7 +50,7 @@ void ElecState::calculate_weights(void)
4950 GlobalC::kv.isk );
5051 Occupy::iweights (GlobalC::kv.nks ,
5152 GlobalC::kv.wk ,
52- GlobalV::NBANDS ,
53+ nbands ,
5354 GlobalC::ucell.magnet .get_neldw (),
5455 ekb_tmp,
5556 GlobalC::en.ef_dw ,
@@ -63,7 +64,7 @@ void ElecState::calculate_weights(void)
6364 // -1 means don't need to consider spin.
6465 Occupy::iweights (GlobalC::kv.nks ,
6566 GlobalC::kv.wk ,
66- GlobalV::NBANDS ,
67+ nbands ,
6768 GlobalC::CHR.nelec ,
6869 ekb_tmp,
6970 this ->ef ,
@@ -77,7 +78,7 @@ void ElecState::calculate_weights(void)
7778 ModuleBase::WARNING_QUIT (" calculate_weights" , " not implemented yet,coming soon!" );
7879 // if(my_rank == 0)
7980 // {
80- // tweights(GlobalC::kv.nkstot, nspin, GlobalV::NBANDS , GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
81+ // tweights(GlobalC::kv.nkstot, nspin, nbands , GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
8182 // this->ef, this->wg);
8283 // }
8384 }
@@ -89,7 +90,7 @@ void ElecState::calculate_weights(void)
8990 double demet_dw = 0.0 ;
9091 Occupy::gweights (GlobalC::kv.nks ,
9192 GlobalC::kv.wk ,
92- GlobalV::NBANDS ,
93+ nbands ,
9394 GlobalC::ucell.magnet .get_nelup (),
9495 Occupy::gaussian_parameter,
9596 Occupy::gaussian_type,
@@ -101,7 +102,7 @@ void ElecState::calculate_weights(void)
101102 GlobalC::kv.isk );
102103 Occupy::gweights (GlobalC::kv.nks ,
103104 GlobalC::kv.wk ,
104- GlobalV::NBANDS ,
105+ nbands ,
105106 GlobalC::ucell.magnet .get_neldw (),
106107 Occupy::gaussian_parameter,
107108 Occupy::gaussian_type,
@@ -118,7 +119,7 @@ void ElecState::calculate_weights(void)
118119 // -1 means is no related to spin.
119120 Occupy::gweights (GlobalC::kv.nks ,
120121 GlobalC::kv.wk ,
121- GlobalV::NBANDS ,
122+ nbands ,
122123 GlobalC::CHR.nelec ,
123124 Occupy::gaussian_parameter,
124125 Occupy::gaussian_type,
@@ -140,7 +141,7 @@ void ElecState::calculate_weights(void)
140141 this ->ef = -1.0e+20 ;
141142 for (int ik = 0 ; ik < GlobalC::kv.nks ; ik++)
142143 {
143- for (int ibnd = 0 ; ibnd < GlobalV::NBANDS ; ibnd++)
144+ for (int ibnd = 0 ; ibnd < nbands ; ibnd++)
144145 {
145146 if (this ->wg (ik, ibnd) > 0.0 )
146147 {
@@ -161,7 +162,7 @@ void ElecState::calculate_weights(void)
161162 double etop = ekb_tmp[0 ][0 ];
162163 for (int ik = 0 ; ik < GlobalC::kv.nks ; ik++)
163164 {
164- for (int ib = 0 ; ib < GlobalV::NBANDS ; ib++)
165+ for (int ib = 0 ; ib < nbands ; ib++)
165166 {
166167 ebotom = min (ebotom, ekb_tmp[ik][ib]);
167168 etop = max (etop, ekb_tmp[ik][ib]);
@@ -185,4 +186,65 @@ void ElecState::calculate_weights(void)
185186 return ;
186187}
187188
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+
188250} // namespace elecstate
0 commit comments