ChargeDensity.C 8.83 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21 22 23
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////

#include "ChargeDensity.h"
#include "Basis.h"
#include "Wavefunction.h"
#include "FourierTransform.h"
#include "SlaterDet.h"
24
#include "blas.h" // dasum
Francois Gygi committed
25 26

#include <iomanip>
Francois Gygi committed
27
#include <algorithm> // fill
28
#include <functional>
29
#include <fstream>
Francois Gygi committed
30 31 32 33
using namespace std;

////////////////////////////////////////////////////////////////////////////////
ChargeDensity::ChargeDensity(const Wavefunction& wf) : ctxt_(wf.context()),
34
wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
Francois Gygi committed
35
{
36
  vbasis_ = new Basis(vcomm_, D3vector(0,0,0));
37
  vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
38
  const Basis& vb = *vbasis_;
39

Francois Gygi committed
40
  // define vft_, FT on vbasis context for transforming the density
41

42 43 44 45 46
  // 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;
Francois Gygi committed
47 48 49
  while (!vbasis_->factorizable(np0v)) np0v += 2;
  while (!vbasis_->factorizable(np1v)) np1v += 2;
  while (!vbasis_->factorizable(np2v)) np2v += 2;
Francois Gygi committed
50
#if 0
51 52 53 54 55 56 57
  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
Francois Gygi committed
58
  vft_ = new FourierTransform(*vbasis_,np0v,np1v,np2v);
59

60
  total_charge_.resize(wf.nspin());
Francois Gygi committed
61 62 63 64 65 66 67 68
  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());
69

Francois Gygi committed
70 71 72 73 74 75
  // 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++ )
    {
76
      const Basis& wb = wf.sd(ispin,ikp)->basis();
Francois Gygi committed
77
#if DEBUG
78 79 80 81 82
      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;
Francois Gygi committed
83
#endif
84 85 86 87 88 89 90 91 92
      // 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);
Francois Gygi committed
93 94
    }
  }
95 96
  // initialize core density ptr to null ptr
  rhocore_r = 0;
Francois Gygi committed
97 98 99 100 101
}

////////////////////////////////////////////////////////////////////////////////
ChargeDensity::~ChargeDensity(void)
{
Francois Gygi committed
102 103
  delete vbasis_;
  delete vft_;
Francois Gygi committed
104
  for ( int ikp = 0; ikp < ft_.size(); ikp++ )
Francois Gygi committed
105 106 107 108 109 110 111 112 113 114
    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 )
    {
115 116
      string s = "name=\"" + (*i).first + "\"";
      cout << "<timing " << left << setw(22) << s
117 118
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
119
           << endl;
Francois Gygi committed
120 121
    }
  }
Francois Gygi committed
122
}
Francois Gygi committed
123

Francois Gygi committed
124 125 126 127 128
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_density(void)
{
  assert(rhor.size() == wf_.nspin());
  const double omega = vbasis_->cell().volume();
129

Francois Gygi committed
130 131
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
Francois Gygi committed
132 133 134 135
    assert(rhor[ispin].size() == vft_->np012loc() );
    assert(rhotmp.size() == vft_->np012loc() );
    const int rhor_size = rhor[ispin].size();
    tmap["charge_compute"].start();
Francois Gygi committed
136
    fill(rhor[ispin].begin(),rhor[ispin].end(),0.0);
Francois Gygi committed
137
    for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
Francois Gygi committed
138
    {
139
      assert(rhor[ispin].size()==ft_[ikp]->np012loc());
Francois Gygi committed
140 141
      wf_.sd(ispin,ikp)->compute_density(*ft_[ikp],
          wf_.weight(ikp), &rhor[ispin][0]);
Francois Gygi committed
142
    }
Francois Gygi committed
143 144 145 146 147 148 149 150 151 152 153 154 155
    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();
Francois Gygi committed
156
    #pragma omp parallel for reduction(+:sum)
Francois Gygi committed
157 158 159 160 161 162 163 164 165 166 167
    for ( int i = 0; i < rhor_size; i++ )
    {
      const double prh = prhor[i];
      sum += prh;
      rhotmp[i] = complex<double>(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();
168
    total_charge_[ispin] = sum;
Francois Gygi committed
169 170 171 172

    tmap["charge_vft"].start();
    vft_->forward(&rhotmp[0],&rhog[ispin][0]);
    tmap["charge_vft"].stop();
173 174 175 176 177 178

    // add core correction charge
    if ( rhocore_r )
      for ( int i = 0; i < rhor_size; i++ )
        rhor[ispin][i] += rhocore_r[i];

Francois Gygi committed
179 180
  }
}
181

Francois Gygi committed
182 183 184 185 186 187 188 189
////////////////////////////////////////////////////////////////////////////////
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;
190

Francois Gygi committed
191 192
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
Francois Gygi committed
193 194 195 196 197 198 199
    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];
200 201 202 203 204 205 206
    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
Francois Gygi committed
207
    {
208 209 210
      #pragma omp parallel for
      for ( int i = 0; i < rhor_size; i++ )
        prhor[i] = rhotmp[i].real() * omega_inv;
Francois Gygi committed
211
    }
212 213 214 215 216 217 218 219 220 221 222 223

    // 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;
Francois Gygi committed
224 225
  }
}
He Ma committed
226

227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
////////////////////////////////////////////////////////////////////////////////
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;
}

He Ma committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_rhog(void)
{
  // recalculate rhog from rhor
  assert(rhor.size() == wf_.nspin());
  const double omega = vbasis_->cell().volume();
  assert(omega!=0.0);

  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    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++ )
        rhotmp[i] = complex<double> ( omega * (prhor[i] - rhocore_r[i]), 0);
    }
    else
    {
      #pragma omp parallel for
      for ( int i = 0; i < rhor_size; i++ )
        rhotmp[i] = complex<double> ( omega * prhor[i], 0);
    }

    assert(rhotmp.size() == vft_->np012loc() );

    vft_->forward(&rhotmp[0],&rhog[ispin][0]);
  }
284
}