////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or .
//
////////////////////////////////////////////////////////////////////////////////
//
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
#include "ChargeDensity.h"
#include "Basis.h"
#include "Wavefunction.h"
#include "FourierTransform.h"
#include "SlaterDet.h"
#include "blas.h" // dasum
#include
#include // fill
#include
using namespace std;
////////////////////////////////////////////////////////////////////////////////
ChargeDensity::ChargeDensity(const Wavefunction& wf) : ctxt_(wf.context()),
wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
{
vbasis_ = new Basis(vcomm_, D3vector(0,0,0));
vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
const Basis& vb = *vbasis_;
// define vft_, FT on vbasis context for transforming the density
// add 2 to grid size to avoid aliasing when using non-zero k-points
// adding 1 would suffice, but add 2 to keep even numbers
int np0v = vbasis_->np(0)+2;
int np1v = vbasis_->np(1)+2;
int np2v = vbasis_->np(2)+2;
while (!vbasis_->factorizable(np0v)) np0v += 2;
while (!vbasis_->factorizable(np1v)) np1v += 2;
while (!vbasis_->factorizable(np2v)) np2v += 2;
#if 0
cout << " ChargeDensity: vbasis: " << endl;
cout << " idxmin: " << vb.idxmin(0) << "/" << vb.idxmin(1)
<< "/" << vb.idxmin(2) << endl;
cout << " idxmax: " << vb.idxmax(0) << "/" << vb.idxmax(1)
<< "/" << vb.idxmax(2) << endl;
cout << " vft grid: " << np0v << "/" << np1v << "/" << np2v << endl;
#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++ )
{
rhor[ispin].resize(vft_->np012loc());
rhog[ispin].resize(vbasis_->localsize());
}
rhotmp.resize(vft_->np012loc());
// FT for interpolation of wavefunctions on the fine grid
ft_.resize(wf.nkp());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const Basis& wb = wf.sd(ispin,ikp)->basis();
#if DEBUG
cout << " ChargeDensity: ikp=" << ikp << " wbasis: " << endl;
cout << " idxmin: " << wb.idxmin(0) << "/" << wb.idxmin(1)
<< "/" << wb.idxmin(2) << endl;
cout << " idxmax: " << wb.idxmax(0) << "/" << wb.idxmax(1)
<< "/" << wb.idxmax(2) << endl;
#endif
// check that no aliasing error can occur
assert(2*np0v > 2*wb.idxmax(0)+vb.idxmax(0));
assert(2*np1v > 2*wb.idxmax(1)+vb.idxmax(1));
assert(2*np2v > 2*wb.idxmax(2)+vb.idxmax(2));
assert(2*np0v > -2*wb.idxmin(0)-vb.idxmin(0));
assert(2*np1v > -2*wb.idxmin(1)-vb.idxmin(1));
assert(2*np2v > -2*wb.idxmin(2)-vb.idxmin(2));
ft_[ikp] = new FourierTransform(wb,np0v,np1v,np2v);
}
}
// initialize core density ptr to null ptr
rhocore_r = 0;
}
////////////////////////////////////////////////////////////////////////////////
ChargeDensity::~ChargeDensity(void)
{
delete vbasis_;
delete vft_;
for ( int ikp = 0; ikp < ft_.size(); ikp++ )
delete ft_[ikp];
for ( TimerMap::iterator i = tmap.begin(); i != tmap.end(); i++ )
{
double time = (*i).second.real();
double tmin = time;
double tmax = time;
ctxt_.dmin(1,1,&tmin,1);
ctxt_.dmax(1,1,&tmax,1);
if ( ctxt_.myproc()==0 )
{
string s = "name=\"" + (*i).first + "\"";
cout << ""
<< endl;
}
}
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_density(void)
{
assert(rhor.size() == wf_.nspin());
const double omega = vbasis_->cell().volume();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
const int rhor_size = rhor[ispin].size();
tmap["charge_compute"].start();
fill(rhor[ispin].begin(),rhor[ispin].end(),0.0);
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
assert(rhor[ispin].size()==ft_[ikp]->np012loc());
wf_.sd(ispin,ikp)->compute_density(*ft_[ikp],
wf_.weight(ikp), &rhor[ispin][0]);
}
tmap["charge_compute"].stop();
// sum over kpoints: sum along rows of kpcontext
tmap["charge_rowsum"].start();
wf_.kpcontext()->dsum('r',vft_->np012loc(),1,
&rhor[ispin][0],vft_->np012loc());
tmap["charge_rowsum"].stop();
// check integral of charge density
// compute Fourier coefficients of the charge density
double sum = 0.0;
const double *const prhor = &rhor[ispin][0];
tmap["charge_integral"].start();
#pragma omp parallel for reduction(+:sum)
for ( int i = 0; i < rhor_size; i++ )
{
const double prh = prhor[i];
sum += prh;
rhotmp[i] = complex(omega * prh, 0.0);
}
sum *= omega / vft_->np012();
// sum on all indices except spin: sum along columns of spincontext
wf_.spincontext()->dsum('c',1,1,&sum,1);
tmap["charge_integral"].stop();
total_charge_[ispin] = sum;
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
tmap["charge_vft"].stop();
// add core correction charge
if ( rhocore_r )
for ( int i = 0; i < rhor_size; i++ )
rhor[ispin][i] += rhocore_r[i];
}
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_rhor(void)
{
// recalculate rhor from rhog
assert(rhor.size() == wf_.nspin());
const double omega = vbasis_->cell().volume();
assert(omega!=0.0);
const double omega_inv = 1.0 / omega;
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
vft_->backward(&rhog[ispin][0],&rhotmp[0]);
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
if ( rhocore_r )
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = ( rhotmp[i].real() + rhocore_r[i] ) * omega_inv;
}
else
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = rhotmp[i].real() * omega_inv;
}
// integral of the charge density
tmap["charge_integral"].start();
int ione=1;
int n = rhor_size;
double sum = dasum(&n,prhor,&ione);
sum *= omega / vft_->np012();
// sum on all indices except spin: sum along columns of spincontext
wf_.spincontext()->dsum('c',1,1,&sum,1);
tmap["charge_integral"].stop();
total_charge_[ispin] = sum;
}
}
////////////////////////////////////////////////////////////////////////////////
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 << " "
<< setprecision(8) << total_charge(ispin)
<< " \n";
}
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const ChargeDensity& cd )
{
cd.print(os);
return os;
}