1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // ChargeDensity.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "MPIdata.h"
20 #include "ChargeDensity.h"
21 #include "Basis.h"
22 #include "Wavefunction.h"
23 #include "FourierTransform.h"
24 #include "SlaterDet.h"
25 #include "blas.h" // dasum
26
27 #include <iomanip>
28 #include <algorithm> // fill
29 #include <functional>
30 #include <fstream>
31 using namespace std;
32
33 ////////////////////////////////////////////////////////////////////////////////
ChargeDensity(const Wavefunction & wf)34 ChargeDensity::ChargeDensity(const Wavefunction& wf) : wf_(wf)
35 {
36 vbasis_ = new Basis(MPIdata::g_comm(), D3vector(0,0,0));
37 vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
38 const Basis& vb = *vbasis_;
39
40 // define vft_, FT on vbasis context for transforming the density
41
42 // add 2 to grid size to avoid aliasing when using non-zero k-points
43 // adding 1 would suffice, but add 2 to keep even numbers
44 int np0v = vb.np(0)+2;
45 int np1v = vb.np(1)+2;
46 int np2v = vb.np(2)+2;
47 while (!vb.factorizable(np0v)) np0v += 2;
48 while (!vb.factorizable(np1v)) np1v += 2;
49 while (!vb.factorizable(np2v)) np2v += 2;
50 #ifdef DEBUG
51 cout << MPIdata::rank() << ": ChargeDensity: vbasis: " << endl;
52 cout << MPIdata::rank() << ": idxmin: "
53 << vb.idxmin(0) << "/" << vb.idxmin(1) << "/" << vb.idxmin(2) << endl;
54 cout << MPIdata::rank() << ": idxmax: "
55 << vb.idxmax(0) << "/" << vb.idxmax(1) << "/" << vb.idxmax(2) << endl;
56 cout << MPIdata::rank() << ": vft grid: "
57 << np0v << "/" << np1v << "/" << np2v << endl;
58 #endif
59 vft_ = new FourierTransform(vb,np0v,np1v,np2v);
60 total_charge_.resize(wf.nspin());
61 rhor.resize(wf.nspin());
62 rhog.resize(wf.nspin());
63 for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
64 {
65 rhor[ispin].resize(vft_->np012loc());
66 rhog[ispin].resize(vb.localsize());
67 }
68 rhotmp.resize(vft_->np012loc());
69
70 #ifdef DEBUG
71 cout << MPIdata::rank() << ": ChargeDensity::ctor: wf.nsp_loc()="
72 << wf.nsp_loc() << " wf.nkp_loc()=" << wf.nkp_loc() << endl;
73 #endif
74
75 // FT for interpolation of wavefunctions on the fine grid
76 if ( wf.nsp_loc() > 0 )
77 {
78 for ( int ikp_loc = 0; ikp_loc < wf.nkp_loc(); ++ikp_loc )
79 {
80 // use basis of isp_loc==0
81 const Basis& wb = wf.sd(0,ikp_loc)->basis();
82 #ifdef DEBUG
83 cout << MPIdata::rank() << ": ChargeDensity::ctor: ikp_loc="
84 << ikp_loc << " wbasis: " << endl;
85 cout << MPIdata::rank() << ": idxmin: " << wb.idxmin(0)
86 << "/" << wb.idxmin(1) << "/" << wb.idxmin(2) << endl;
87 cout << MPIdata::rank() << ": idxmax: " << wb.idxmax(0)
88 << "/" << wb.idxmax(1) << "/" << wb.idxmax(2) << endl;
89 #endif
90 // check that no aliasing error can occur
91 assert(2*np0v > 2*wb.idxmax(0)+vb.idxmax(0));
92 assert(2*np1v > 2*wb.idxmax(1)+vb.idxmax(1));
93 assert(2*np2v > 2*wb.idxmax(2)+vb.idxmax(2));
94 assert(2*np0v > -2*wb.idxmin(0)-vb.idxmin(0));
95 assert(2*np1v > -2*wb.idxmin(1)-vb.idxmin(1));
96 assert(2*np2v > -2*wb.idxmin(2)-vb.idxmin(2));
97
98 ft_.push_back(new FourierTransform(wb,np0v,np1v,np2v));
99 }
100 }
101 // perform reset operation on timers to introduce them in tmap
102 // this is necessary to have a tmap entry even if nsp_loc or nkp_loc == 0
103 tmap["charge_compute"].reset();
104 tmap["charge_rowsum"].reset();
105 tmap["charge_integral"].reset();
106 tmap["charge_vft"].reset();
107 tmap["update_taur"].reset();
108 }
109
110 ////////////////////////////////////////////////////////////////////////////////
~ChargeDensity(void)111 ChargeDensity::~ChargeDensity(void)
112 {
113 delete vbasis_;
114 delete vft_;
115 for ( int ikp_loc = 0; ikp_loc < ft_.size(); ++ikp_loc )
116 delete ft_[ikp_loc];
117
118 for ( TimerMap::iterator i = tmap.begin(); i != tmap.end(); i++ )
119 {
120 double time = i->second.real();
121 double tmin, tmax;
122 MPI_Reduce(&time,&tmin,1,MPI_DOUBLE,MPI_MIN,0,MPIdata::comm());
123 MPI_Reduce(&time,&tmax,1,MPI_DOUBLE,MPI_MAX,0,MPIdata::comm());
124 if ( MPIdata::onpe0() && (tmax > 0.0) )
125 {
126 string s = "name=\"" + (*i).first + "\"";
127 cout << "<timing " << left << setw(22) << s
128 << " min=\"" << setprecision(3) << tmin << "\""
129 << " max=\"" << setprecision(3) << tmax << "\"/>"
130 << endl;
131 }
132 }
133 }
134
135 ////////////////////////////////////////////////////////////////////////////////
update_density(void)136 void ChargeDensity::update_density(void)
137 {
138 assert(rhor.size() == wf_.nspin());
139 for ( int ispin = 0; ispin < wf_.nspin(); ++ispin )
140 fill(rhor[ispin].begin(),rhor[ispin].end(),0.0);
141
142 for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
143 {
144 const int ispg = wf_.isp_global(isp_loc);
145 assert(rhor[ispg].size() == vft_->np012loc() );
146 assert(rhotmp.size() == vft_->np012loc() );
147
148 tmap["charge_compute"].start();
149 for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
150 {
151 assert(ft_[ikp_loc]);
152 const int ikpg = wf_.ikp_global(ikp_loc);
153 assert(rhor[ispg].size()==ft_[ikp_loc]->np012loc());
154 wf_.sd(isp_loc,ikp_loc)->compute_density(*ft_[ikp_loc],
155 wf_.weight(ikpg), &rhor[ispg][0]);
156 }
157 tmap["charge_compute"].stop();
158 }
159 // rhor now contains local contributions from this task
160
161 for ( int ispin = 0; ispin < wf_.nspin(); ++ispin )
162 {
163 // sum over kpoints, states and spins
164 tmap["charge_rowsum"].start();
165 vector<double> tmpr(vft_->np012loc());
166 MPI_Allreduce(&rhor[ispin][0],&tmpr[0],vft_->np012loc(),
167 MPI_DOUBLE,MPI_SUM,MPIdata::st_kp_sp_comm());
168 rhor[ispin] = tmpr;
169 tmap["charge_rowsum"].stop();
170
171 // check integral of charge density
172 // compute Fourier coefficients of the charge density
173 double sum = 0.0;
174 const double *const prhor = &rhor[ispin][0];
175 tmap["charge_integral"].start();
176 const int rhor_size = rhor[ispin].size();
177 const double omega = vbasis_->cell().volume();
178 #pragma omp parallel for reduction(+:sum)
179 for ( int i = 0; i < rhor_size; i++ )
180 {
181 const double prh = prhor[i];
182 sum += prh;
183 rhotmp[i] = complex<double>(omega * prh, 0.0);
184 }
185 sum *= omega / vft_->np012();
186 double tsum = 0.0;
187 // sum over g_comm and sp_comm
188 MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::g_comm());
189 sum = tsum;
190 tmap["charge_integral"].stop();
191 total_charge_[ispin] = sum;
192
193 // compute rhog from rhotmp
194 tmap["charge_vft"].start();
195 vft_->forward(&rhotmp[0],&rhog[ispin][0]);
196 tmap["charge_vft"].stop();
197 }
198 }
199
200 ////////////////////////////////////////////////////////////////////////////////
update_rhor(void)201 void ChargeDensity::update_rhor(void)
202 {
203 // recalculate rhor from rhog
204 assert(rhor.size() == wf_.nspin());
205 const double omega = vbasis_->cell().volume();
206 assert(omega!=0.0);
207 const double omega_inv = 1.0 / omega;
208
209 for ( int ispin = 0; ispin < wf_.nspin(); ++ispin )
210 {
211 assert(rhor[ispin].size() == vft_->np012loc() );
212 assert(rhotmp.size() == vft_->np012loc() );
213
214 vft_->backward(&rhog[ispin][0],&rhotmp[0]);
215
216 const int rhor_size = rhor[ispin].size();
217 double *const prhor = &rhor[ispin][0];
218 #pragma omp parallel for
219 for ( int i = 0; i < rhor_size; i++ )
220 prhor[i] = rhotmp[i].real() * omega_inv;
221
222 // integral of the charge density
223 tmap["charge_integral"].start();
224 int ione=1;
225 int n = rhor_size;
226 double sum = dasum(&n,prhor,&ione);
227 sum *= omega / vft_->np012();
228
229 // sum over g_comm
230 double tsum = 0.0;
231 MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::g_comm());
232 tmap["charge_integral"].stop();
233 total_charge_[ispin] = tsum;
234 }
235 }
236
237 ////////////////////////////////////////////////////////////////////////////////
update_taur(double * taur) const238 void ChargeDensity::update_taur(double* taur) const
239 {
240 vector<double> tmpr(vft_->np012loc(),0.0);
241 tmap["update_taur"].start();
242 for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
243 {
244 for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
245 {
246 assert(ft_[ikp_loc]);
247 const int ikpg = wf_.ikp_global(ikp_loc);
248 wf_.sd(isp_loc,ikp_loc)->compute_tau(*ft_[ikp_loc],
249 wf_.weight(ikpg), &tmpr[0]);
250 }
251 }
252 MPI_Allreduce(&tmpr[0],taur,vft_->np012loc(),
253 MPI_DOUBLE,MPI_SUM,MPIdata::st_kp_sp_comm());
254 tmap["update_taur"].stop();
255
256 // stop if computing taur with NLCCs
257 if ( !rhocore_r.empty() )
258 assert(!"ChargeDensity: Cannot compute taur with NLCCs");
259 }
260
261 ////////////////////////////////////////////////////////////////////////////////
update_taur(double * taur_up,double * taur_dn) const262 void ChargeDensity::update_taur(double* taur_up, double* taur_dn) const
263 {
264 assert(wf_.nspin()==2);
265 vector<double> tmpr(vft_->np012loc());
266 tmap["update_taur"].start();
267 for ( int ispin = 0; ispin < wf_.nspin(); ++ispin )
268 {
269 memset( (void*)&tmpr[0], 0, vft_->np012loc()*sizeof(double) );
270 const int isp_loc = wf_.isp_local(ispin);
271 if ( isp_loc >= 0 )
272 {
273 for ( int ikp_loc = 0; ikp_loc < wf_.nkp_loc(); ++ikp_loc )
274 {
275 assert(ft_[ikp_loc]);
276 const int ikpg = wf_.ikp_global(ikp_loc);
277 wf_.sd(isp_loc,ikp_loc)->compute_tau(*ft_[ikp_loc],
278 wf_.weight(ikpg), &tmpr[0]);
279 }
280 }
281 // tmpr now contains local taur contrib
282 // reduce to have both taur_up and taur_dn on all tasks
283 if ( ispin == 0 )
284 {
285 MPI_Allreduce(&tmpr[0],&taur_up[0],vft_->np012loc(),
286 MPI_DOUBLE,MPI_SUM,MPIdata::st_kp_sp_comm());
287 }
288 else
289 {
290 MPI_Allreduce(&tmpr[0],&taur_dn[0],vft_->np012loc(),
291 MPI_DOUBLE,MPI_SUM,MPIdata::st_kp_sp_comm());
292 }
293 }
294 tmap["update_taur"].stop();
295
296 // stop if computing taur with NLCCs
297 if ( !rhocore_r.empty() )
298 assert(!"ChargeDensity: Cannot compute taur with NLCCs");
299 }
300
301 ////////////////////////////////////////////////////////////////////////////////
total_charge(void) const302 double ChargeDensity::total_charge(void) const
303 {
304 assert((wf_.nspin()==1)||(wf_.nspin()==2));
305 if ( wf_.nspin() == 1 )
306 return total_charge_[0];
307 else
308 return total_charge_[0] + total_charge_[1];
309 }
310
311 ////////////////////////////////////////////////////////////////////////////////
print(ostream & os) const312 void ChargeDensity::print(ostream& os) const
313 {
314 os.setf(ios::fixed,ios::floatfield);
315 os.setf(ios::right,ios::adjustfield);
316 for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
317 os << " <electronic_charge ispin=\"" << ispin << "\"> "
318 << setprecision(8) << total_charge(ispin)
319 << " </electronic_charge>\n";
320 }
321
322 ////////////////////////////////////////////////////////////////////////////////
operator <<(ostream & os,const ChargeDensity & cd)323 ostream& operator<< ( ostream& os, const ChargeDensity& cd )
324 {
325 cd.print(os);
326 return os;
327 }
328
329 ////////////////////////////////////////////////////////////////////////////////
update_rhog(void)330 void ChargeDensity::update_rhog(void)
331 {
332 // recalculate rhog from rhor
333 assert(rhor.size() == wf_.nspin());
334 const double omega = vbasis_->cell().volume();
335 assert(omega!=0.0);
336
337 for ( int isp_loc = 0; isp_loc < wf_.nsp_loc(); ++isp_loc )
338 {
339 const int ispg = wf_.isp_global(isp_loc);
340 const int rhor_size = rhor[ispg].size();
341 double *const prhor = &rhor[ispg][0];
342 #pragma omp parallel for
343 for ( int i = 0; i < rhor_size; i++ )
344 rhotmp[i] = complex<double> ( omega * prhor[i], 0);
345
346 assert(rhotmp.size() == vft_->np012loc() );
347
348 vft_->forward(&rhotmp[0],&rhog[ispg][0]);
349 }
350 }
351