Commit 26c34952 by Francois Gygi

Update printing of electronic charge & net charge

parent c102a324
......@@ -346,6 +346,9 @@ void BOSampleStepper::step(int niter)
wkerker[i] = 1.0;
}
if ( onpe0 )
cout << "<net_charge> " << atoms.nel()-wf.nel() << " </net_charge>\n";
// Next line: special case of niter=0: compute GS only
for ( int iter = 0; iter < max(niter,1); iter++ )
{
......@@ -376,6 +379,7 @@ void BOSampleStepper::step(int niter)
if ( onpe0 )
{
cout << cd_;
cout << ef_;
if ( ef_.el_enth() )
cout << *ef_.el_enth();
......@@ -756,6 +760,9 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
if ( onpe0 )
cout << cd_;
// charge mixing
if ( nite_ > 0 )
{
......@@ -1173,10 +1180,10 @@ void BOSampleStepper::step(int niter)
tmap["energy"].stop();
if ( onpe0 )
{
cout << cd_;
cout << ef_;
if ( ef_.el_enth() )
cout << *ef_.el_enth();
}
}
......
......@@ -55,6 +55,7 @@ wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
#endif
vft_ = new FourierTransform(*vbasis_,np0v,np1v,np2v);
total_charge_.resize(wf.nspin());
rhor.resize(wf.nspin());
rhog.resize(wf.nspin());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
......@@ -162,13 +163,7 @@ void ChargeDensity::update_density(void)
// sum on all indices except spin: sum along columns of spincontext
wf_.spincontext()->dsum('c',1,1,&sum,1);
tmap["charge_integral"].stop();
if ( ctxt_.onpe0() )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " total_electronic_charge: " << setprecision(8) << sum
<< endl;
}
total_charge_[ispin] = sum;
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
......@@ -214,3 +209,31 @@ void ChargeDensity::update_rhor(void)
}
}
}
////////////////////////////////////////////////////////////////////////////////
double ChargeDensity::total_charge(void) const
{
assert((wf_.nspin()==1)||(wf_.nspin()==2));
if ( wf_.nspin() == 1 )
return total_charge_[0];
else
return total_charge_[0] + total_charge_[1];
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::print(ostream& os) const
{
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
os << " <electronic_charge ispin=\"" << ispin << "\"> "
<< setprecision(8) << total_charge(ispin)
<< " </electronic_charge>\n";
}
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const ChargeDensity& cd )
{
cd.print(os);
return os;
}
......@@ -44,6 +44,7 @@ class ChargeDensity
FourierTransform* vft_;
std::vector<FourierTransform*> ft_; // ft_[ikp];
std::valarray<std::complex<double> > rhotmp;
std::vector<double> total_charge_;
public:
......@@ -61,8 +62,13 @@ class ChargeDensity
Basis* vbasis(void) const { return vbasis_; }
FourierTransform* vft(void) const { return vft_; }
FourierTransform* ft(int ikp) const { return ft_[ikp]; }
double total_charge(int ispin) const { return total_charge_[ispin]; }
double total_charge(void) const;
void print(std::ostream& os) const;
ChargeDensity(const Wavefunction& wf);
~ChargeDensity();
};
std::ostream& operator << ( std::ostream& os, const ChargeDensity& cd );
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment