ChargeDensity.C 5.05 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
6
// $Id: ChargeDensity.C,v 1.8 2004-09-14 22:24:11 fgygi Exp $
Francois Gygi committed
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 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 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

#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()),
wf_(wf) 
{
  // determine which k points are stored on the sd context to which 
  // this pe belongs
  int myspin = -1;
  int ikplast = -1;
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
      if ( wf.sd(ispin,ikp) != 0 )
      {
        myspin = ispin;
        ikplast = ikp;
      }
    }
  }
  // the last kpoint on this process is ikplast
  
  // create a vbasis
  // Note: avoid creating more than one vcontext if there are 
  // more than one kpoint per process
  if ( wf.sd(myspin,ikplast) != 0 )
  {
    // Note: Basis defined on columns of vcontext
    vcontext_ = new Context(*wf_.spincontext(myspin),'c',
      wf_.spincontext(myspin)->mycol());
    vbasis_ = new Basis(*vcontext_, D3vector(0,0,0));
    vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
    //cout << vcontext_->mype() << ": vcontext = " << *vcontext_ << endl;
  }
  else
  {
    vcontext_ = 0;
    vbasis_ = 0;
  }
    
  // 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)
{
  delete vcontext_;
  delete vbasis_;
  delete vft_;
  for ( int ikp = 0; ikp < ft_.size(); ikp++ )
    delete ft_[ikp];
}
Francois Gygi committed
96

Francois Gygi committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
////////////////////////////////////////////////////////////////////////////////
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() );
      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]);
 
        }
      }
    
      // sum across rows of spincontext[ispin] context
      //tmap["dsum"].start();
      wf_.spincontext(ispin)->dsum('r',vft_->np012loc(),1,
        &rhor[ispin][0],vft_->np012loc());
      //tmap["dsum"].stop();
      
      // check integral of charge density
Francois Gygi committed
127
      // compute Fourier coefficients of the charge density
Francois Gygi committed
128
      double sum = 0.0;
Francois Gygi committed
129 130 131 132 133 134 135 136 137
      const int rhor_size = rhor[ispin].size();
      const double *const prhor = &rhor[ispin][0];
      #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
138 139 140 141 142 143 144
      sum *= omega / vft_->np012();

      wf_.spincontext(ispin)->dsum('c',1,1,&sum,1);
      if ( ctxt_.onpe0() )
      {
        cout.setf(ios::fixed,ios::floatfield);
        cout.setf(ios::right,ios::adjustfield);
Francois Gygi committed
145 146
        cout << "  <!-- total_electronic_charge: " << setprecision(8) << sum 
             << " -->" << endl;
Francois Gygi committed
147
      }
Francois Gygi committed
148
        
Francois Gygi committed
149 150 151 152
      vft_->forward(&rhotmp[0],&rhog[ispin][0]);
    }
  }
}
Francois Gygi committed
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
////////////////////////////////////////////////////////////////////////////////
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
181 182