ChargeDensity.C 5.05 KB
 Francois Gygi committed Oct 02, 2003 1 2 3 4 5 //////////////////////////////////////////////////////////////////////////////// // // ChargeDensity.C // ////////////////////////////////////////////////////////////////////////////////  Francois Gygi committed Sep 14, 2004 6 // $Id: ChargeDensity.C,v 1.8 2004-09-14 22:24:11 fgygi Exp$  Francois Gygi committed Oct 02, 2003 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 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 Sep 14, 2004 96   Francois Gygi committed Oct 02, 2003 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 Jun 01, 2004 127  // compute Fourier coefficients of the charge density  Francois Gygi committed Oct 02, 2003 128  double sum = 0.0;  Francois Gygi committed Jun 01, 2004 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(omega * prh, 0.0); }  Francois Gygi committed Oct 02, 2003 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 Aug 11, 2004 145 146  cout << " " << endl;  Francois Gygi committed Oct 02, 2003 147  }  Francois Gygi committed Jun 01, 2004 148   Francois Gygi committed Oct 02, 2003 149 150 151 152  vft_->forward(&rhotmp[0],&rhog[ispin][0]); } } }  Francois Gygi committed Sep 14, 2004 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 Oct 02, 2003 181 182