ChargeDensity.C 4.93 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: ChargeDensity.C,v 1.11 2004-11-10 22:30:53 fgygi Exp $
Francois Gygi committed
7 8 9 10 11 12 13 14 15 16 17 18

#include "ChargeDensity.h"
#include "Basis.h"
#include "Wavefunction.h"
#include "FourierTransform.h"
#include "SlaterDet.h"

#include <iomanip>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
ChargeDensity::ChargeDensity(const Wavefunction& wf) : ctxt_(wf.context()),
19 20
wf_(wf), vcontext_(wf.sd(0,0)->basis().context()) 

Francois Gygi committed
21
{
22 23
  vbasis_ = new Basis(vcontext_, D3vector(0,0,0));
  vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
Francois Gygi committed
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
    
  // define vft_, FT on vbasis context for transforming the density
  
  int np0v = vbasis_->np(0);
  int np1v = vbasis_->np(1);
  int np2v = vbasis_->np(2);
  vft_ = new FourierTransform(*vbasis_,np0v,np1v,np2v);
  
  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++ )
    {
      if ( wf.sd(ispin,ikp) != 0 )
        ft_[ikp] =
          new FourierTransform(wf.sd(ispin,ikp)->basis(),np0v,np1v,np2v);
      else
        ft_[ikp] = 0;
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
ChargeDensity::~ChargeDensity(void)
{
Francois Gygi committed
59 60
  delete vbasis_;
  delete vft_;
Francois Gygi committed
61
  for ( int ikp = 0; ikp < ft_.size(); ikp++ )
Francois Gygi committed
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
    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 )
    {
      cout << "<!-- timing "
           << setw(15) << (*i).first
           << " : " << setprecision(3) << setw(9) << tmin
           << " "   << setprecision(3) << setw(9) << tmax << " -->" << endl;
    }
  }
Francois Gygi committed
78
}
Francois Gygi committed
79

Francois Gygi committed
80 81 82 83 84 85 86 87 88 89 90 91
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_density(void)
{
  assert(rhor.size() == wf_.nspin());
  const double omega = vbasis_->cell().volume();
  
  for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
  {
    if ( wf_.spincontext(ispin)->active())
    {
      assert(rhor[ispin].size() == vft_->np012loc() );
      assert(rhotmp.size() == vft_->np012loc() );
Francois Gygi committed
92
      tmap["charge_compute"].start();
Francois Gygi committed
93 94 95 96 97 98 99 100 101 102
      for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
      {
        if ( wf_.sd(ispin,ikp) != 0 && wf_.sdcontext(ispin,ikp)->active() )
        {
      
          wf_.sd(ispin,ikp)->compute_density(*ft_[ikp],
            wf_.weight(ikp), &rhor[ispin][0]);
 
        }
      }
Francois Gygi committed
103
      tmap["charge_compute"].stop();
Francois Gygi committed
104 105
    
      // sum across rows of spincontext[ispin] context
Francois Gygi committed
106
      tmap["charge_rowsum"].start();
Francois Gygi committed
107 108
      wf_.spincontext(ispin)->dsum('r',vft_->np012loc(),1,
        &rhor[ispin][0],vft_->np012loc());
Francois Gygi committed
109
      tmap["charge_rowsum"].stop();
Francois Gygi committed
110 111
      
      // check integral of charge density
Francois Gygi committed
112
      // compute Fourier coefficients of the charge density
Francois Gygi committed
113
      double sum = 0.0;
Francois Gygi committed
114 115
      const int rhor_size = rhor[ispin].size();
      const double *const prhor = &rhor[ispin][0];
Francois Gygi committed
116
      tmap["charge_integral"].start();
Francois Gygi committed
117 118 119 120 121 122 123
      #pragma ivdep
      for ( int i = 0; i < rhor_size; i++ )
      {
        const double prh = prhor[i];
        sum += prh;
        rhotmp[i] = complex<double>(omega * prh, 0.0);
      }
Francois Gygi committed
124 125 126
      sum *= omega / vft_->np012();

      wf_.spincontext(ispin)->dsum('c',1,1,&sum,1);
Francois Gygi committed
127
      tmap["charge_integral"].stop();
Francois Gygi committed
128 129 130 131
      if ( ctxt_.onpe0() )
      {
        cout.setf(ios::fixed,ios::floatfield);
        cout.setf(ios::right,ios::adjustfield);
Francois Gygi committed
132 133
        cout << "  <!-- total_electronic_charge: " << setprecision(8) << sum 
             << " -->" << endl;
Francois Gygi committed
134
      }
Francois Gygi committed
135
        
136
      tmap["charge_vft"].start();
Francois Gygi committed
137
      vft_->forward(&rhotmp[0],&rhog[ispin][0]);
138
      tmap["charge_vft"].stop();
Francois Gygi committed
139 140 141
    }
  }
}
Francois Gygi committed
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
////////////////////////////////////////////////////////////////////////////////
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++ )
  {
    if ( wf_.spincontext(ispin)->active())
    {
      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];
      #pragma ivdep
      for ( int i = 0; i < rhor_size; i++ )
      {
        prhor[i] = rhotmp[i].real() * omega_inv;
      }
    }
  }
}
Francois Gygi committed
170 171