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 // Wavefunction.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "MPIdata.h"
20 #include "Wavefunction.h"
21 #include "SlaterDet.h"
22 #include "FourierTransform.h"
23 #include "jacobi.h"
24 #include "SharedFilePtr.h"
25 #include "cout0.h"
26 #include <vector>
27 #include <iomanip>
28 #include <sstream>
29 using namespace std;
30
31 ////////////////////////////////////////////////////////////////////////////////
Wavefunction(const Context & sd_ctxt)32 Wavefunction::Wavefunction(const Context& sd_ctxt) : sd_ctxt_(sd_ctxt), nel_(0),
33 nempty_(0), nspin_(1), deltaspin_(0), ecut_(0.0)
34 {
35 // create a default wavefunction: one k point, k=0
36 kpoint_.resize(1);
37 kpoint_[0] = D3vector(0,0,0);
38 weight_.resize(1);
39 weight_[0] = 1.0;
40 compute_nst();
41 allocate();
42 }
43
44 ////////////////////////////////////////////////////////////////////////////////
Wavefunction(const Wavefunction & wf)45 Wavefunction::Wavefunction(const Wavefunction& wf) : sd_ctxt_(wf.sd_ctxt_),
46 nel_(wf.nel_), nempty_(wf.nempty_), nspin_(wf.nspin_), nst_(wf.nst_),
47 deltaspin_(wf.deltaspin_), cell_(wf.cell_), refcell_(wf.refcell_),
48 ecut_(wf.ecut_), weight_(wf.weight_), kpoint_(wf.kpoint_)
49 {
50 allocate();
51 resize();
52 init();
53 }
54
55 ////////////////////////////////////////////////////////////////////////////////
~Wavefunction()56 Wavefunction::~Wavefunction()
57 {
58 deallocate();
59 }
60
61 ////////////////////////////////////////////////////////////////////////////////
allocate(void)62 void Wavefunction::allocate(void)
63 {
64 // compute local number of kpoints nkp_loc_[ikpb]
65 nkp_loc_.resize(MPIdata::nkpb());
66 const int nkp = kpoint_.size();
67 //cout << MPIdata::rank() << ": nkp=" << nkp << endl;
68 //cout << MPIdata::rank() << ": MPIdata::nkpb()=" << MPIdata::nkpb() << endl;
69 for ( int ikpb = 0; ikpb < MPIdata::nkpb(); ++ikpb )
70 {
71 nkp_loc_[ikpb] = nkp / MPIdata::nkpb() +
72 (ikpb < (nkp % MPIdata::nkpb()) ? 1 : 0);
73 //cout << MPIdata::rank() << ": nkp_loc_[" << ikpb << "]="
74 // << nkp_loc_[ikpb] << endl;
75 }
76
77 // round robin allocation of kpoints
78 ikp_global_.resize(0);
79 for ( int ikpg = MPIdata::ikpb(); ikpg < nkp; ikpg += MPIdata::nkpb() )
80 {
81 ikp_global_.push_back(ikpg);
82 }
83
84 //cout << MPIdata::rank() << ": nkp_loc_[MPIdata::ikpb]="
85 // << nkp_loc_[MPIdata::ikpb()] << endl;
86 //cout << MPIdata::rank() << ": ikp_global_.size()="
87 // << ikp_global_.size() << endl;
88 assert(nkp_loc_[MPIdata::ikpb()] == ikp_global_.size());
89
90 // compute local number of spins nsp_loc_[ispb]
91 nsp_loc_.resize(MPIdata::nspb());
92 for ( int ispb = 0; ispb < MPIdata::nspb(); ++ispb)
93 {
94 nsp_loc_[ispb] = nspin_ / MPIdata::nspb() +
95 (ispb < (nspin_ % MPIdata::nspb()) ? 1 : 0);
96 //cout << MPIdata::rank() << ": nsp_loc_[" << ispb << "]="
97 // << nsp_loc_[ispb] << endl;
98 }
99
100 // round robin allocation of spins
101 isp_global_.resize(0);
102 for ( int ispg = MPIdata::ispb(); ispg < nspin_; ispg += MPIdata::nspb() )
103 {
104 isp_global_.push_back(ispg);
105 }
106
107 assert(nsp_loc_[MPIdata::ispb()] == isp_global_.size());
108
109 // create local SlaterDets
110 sd_.resize(nsp_loc_[MPIdata::ispb()]);
111 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
112 {
113 sd_[isp_loc].resize(nkp_loc_[MPIdata::ikpb()]);
114 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
115 {
116 int ikp = ikp_global_[ikp_loc];
117 sd_[isp_loc][ikp_loc] = new SlaterDet(sd_ctxt_,kpoint_[ikp]);
118 }
119 }
120 }
121
122 ////////////////////////////////////////////////////////////////////////////////
deallocate(void)123 void Wavefunction::deallocate(void)
124 {
125 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
126 {
127 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
128 {
129 delete sd_[isp_loc][ikp_loc];
130 }
131 sd_[isp_loc].resize(0);
132 }
133 }
134
135 ////////////////////////////////////////////////////////////////////////////////
nkp(void) const136 int Wavefunction::nkp(void) const { return kpoint_.size(); }
137
138 ////////////////////////////////////////////////////////////////////////////////
nkp_loc() const139 int Wavefunction::nkp_loc() const { return nkp_loc_[MPIdata::ikpb()]; }
140
141 ////////////////////////////////////////////////////////////////////////////////
ikp_global(int ikp) const142 int Wavefunction::ikp_global(int ikp) const { return ikp_global_[ikp]; }
143
144 ////////////////////////////////////////////////////////////////////////////////
ikp_local(int ikpg) const145 int Wavefunction::ikp_local(int ikpg) const
146 {
147 // return local index of ikpg if hosted on this task, -1 otherwise
148 if ( ( ikpg % MPIdata::nkpb() ) == MPIdata::ikpb() )
149 return ikpg / MPIdata::nkpb();
150 else
151 return -1;
152 }
153
154 ////////////////////////////////////////////////////////////////////////////////
nel() const155 int Wavefunction::nel() const { return nel_; } // total number of electrons
156
157 ////////////////////////////////////////////////////////////////////////////////
nst() const158 int Wavefunction::nst() const
159 {
160 if ( nspin_ == 1 )
161 return nst_[0];
162 else
163 return nst_[0]+nst_[1];
164 }
165
166 ////////////////////////////////////////////////////////////////////////////////
nst(int ispin) const167 int Wavefunction::nst(int ispin) const
168 {
169 assert(ispin >= 0 && ispin < 2);
170 return nst_[ispin];
171 }
172
173 ////////////////////////////////////////////////////////////////////////////////
nempty() const174 int Wavefunction::nempty() const { return nempty_; }
175
176 ////////////////////////////////////////////////////////////////////////////////
nspin() const177 int Wavefunction::nspin() const { return nspin_; }
178
179 ////////////////////////////////////////////////////////////////////////////////
nsp_loc() const180 int Wavefunction::nsp_loc() const { return nsp_loc_[MPIdata::ispb()]; }
181
182 ////////////////////////////////////////////////////////////////////////////////
isp_global(int isp) const183 int Wavefunction::isp_global(int isp) const { return isp_global_[isp]; }
184
185 ////////////////////////////////////////////////////////////////////////////////
isp_local(int ispg) const186 int Wavefunction::isp_local(int ispg) const
187 {
188 // return local index of ispg if hosted on this task, -1 otherwise
189 if ( ( ispg % MPIdata::nspb() ) == MPIdata::ispb() )
190 return ispg / MPIdata::nspb();
191 else
192 return -1;
193 }
194
195 ////////////////////////////////////////////////////////////////////////////////
deltaspin() const196 int Wavefunction::deltaspin() const { return deltaspin_; }
197
198 ////////////////////////////////////////////////////////////////////////////////
entropy(void) const199 double Wavefunction::entropy(void) const
200 {
201 double sum = 0.0;
202 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
203 {
204 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
205 {
206 int ikp = ikp_global_[ikp_loc];
207 sum += weight_[ikp] * sd(isp_loc,ikp_loc)->entropy(nspin_);
208 }
209 }
210 double tsum;
211 MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::kp_sp_comm());
212 return tsum;
213 }
214
215 ////////////////////////////////////////////////////////////////////////////////
resize(void)216 void Wavefunction::resize(void)
217 {
218 try
219 {
220 // resize all SlaterDets using current cell_, refcell_, ecut_, nst_[ispin]
221 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
222 {
223 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
224 {
225 int isp = isp_global_[isp_loc];
226 sd_[isp_loc][ikp_loc]->resize(cell_,refcell_,ecut_,nst_[isp]);
227 }
228 }
229 }
230 catch ( const SlaterDetException& sdex )
231 {
232 cout << " Wavefunction: SlaterDetException during resize: " << endl
233 << sdex.msg << endl;
234 // no resize took place
235 return;
236 }
237 catch ( bad_alloc )
238 {
239 cout << " Wavefunction: insufficient memory for resize operation" << endl;
240 return;
241 }
242 }
243
244 ////////////////////////////////////////////////////////////////////////////////
resize(const UnitCell & cell,const UnitCell & refcell,double ecut)245 void Wavefunction::resize(const UnitCell& cell, const UnitCell& refcell,
246 double ecut)
247 {
248 try
249 {
250 // resize all SlaterDets using cell, refcell, ecut and nst_[ispin]
251 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
252 {
253 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
254 {
255 int isp = isp_global_[isp_loc];
256 sd_[isp_loc][ikp_loc]->resize(cell,refcell,ecut,nst_[isp]);
257 }
258 }
259 cell_ = cell;
260 refcell_ = refcell;
261 ecut_ = ecut;
262 }
263 catch ( const SlaterDetException& sdex )
264 {
265 cout << " Wavefunction: SlaterDetException during resize: " << endl
266 << sdex.msg << endl;
267 // no resize took place
268 return;
269 }
270 catch ( bad_alloc )
271 {
272 cout << " Wavefunction: insufficient memory for resize operation" << endl;
273 return;
274 }
275 }
276
277 ////////////////////////////////////////////////////////////////////////////////
init(void)278 void Wavefunction::init(void)
279 {
280 // initialize all SlaterDets with lowest energy plane waves
281 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
282 {
283 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
284 {
285 sd_[isp_loc][ikp_loc]->init();
286 }
287 }
288 }
289
290 ////////////////////////////////////////////////////////////////////////////////
clear(void)291 void Wavefunction::clear(void)
292 {
293 // initialize all SlaterDets with zero
294 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
295 {
296 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
297 {
298 sd_[isp_loc][ikp_loc]->c().clear();
299 }
300 }
301 }
302
303 ////////////////////////////////////////////////////////////////////////////////
reset(void)304 void Wavefunction::reset(void)
305 {
306 // reset to single kpoint, ecut=0
307 deallocate();
308 nel_ = 0;
309 nempty_ = 0;
310 nspin_ = 1;
311 deltaspin_ = 0;
312 kpoint_.resize(1);
313 kpoint_[0] = D3vector(0,0,0);
314 weight_.resize(1);
315 weight_[0] = 1.0;
316 compute_nst();
317 allocate();
318 }
319
320 ////////////////////////////////////////////////////////////////////////////////
compute_nst(void)321 void Wavefunction::compute_nst(void)
322 {
323 // recompute nst from nel_, deltaspin_, nempty_
324
325 nst_.resize(nspin_);
326 if ( nspin_ == 1 )
327 {
328 nst_[0] = ( nel_ + 1 ) / 2 + nempty_;
329 }
330 else
331 {
332 // nspin == 2
333 nst_[0] = ( nel_ + 1 ) / 2 + deltaspin_ + nempty_;
334 nst_[1] = nel_ / 2 - deltaspin_ + nempty_;
335 }
336 }
337
338 ////////////////////////////////////////////////////////////////////////////////
set_nel(int nel)339 void Wavefunction::set_nel(int nel)
340 {
341 if ( nel == nel_ ) return;
342 if ( nel < 0 )
343 {
344 cout << " Wavefunction::set_nel: nel < 0" << endl;
345 return;
346 }
347
348 nel_ = nel;
349 compute_nst();
350 resize();
351 }
352
353 ////////////////////////////////////////////////////////////////////////////////
set_nempty(int nempty)354 void Wavefunction::set_nempty(int nempty)
355 {
356 if ( nempty == nempty_ ) return;
357 if ( nempty < 0 )
358 {
359 cout << " Wavefunction::set_nempty: negative value" << endl;
360 return;
361 }
362 nempty_ = nempty;
363 compute_nst();
364 update_occ(0.0);
365 resize();
366 }
367
368 ////////////////////////////////////////////////////////////////////////////////
set_nspin(int nspin)369 void Wavefunction::set_nspin(int nspin)
370 {
371 assert(nspin==1 || nspin==2);
372 if ( nspin == nspin_ ) return;
373
374 deallocate();
375 nspin_ = nspin;
376
377 compute_nst();
378 allocate();
379 resize();
380 init();
381 update_occ(0.0);
382 }
383
384 ////////////////////////////////////////////////////////////////////////////////
set_deltaspin(int deltaspin)385 void Wavefunction::set_deltaspin(int deltaspin)
386 {
387 if ( deltaspin == deltaspin_ ) return;
388
389 // check if value of deltaspin would result in nst[1] < 0
390 // nst_[1] = nel_ / 2 - deltaspin_ + nempty_;
391 if ( nel_ / 2 - deltaspin + nempty_ < 0 )
392 {
393 if ( MPIdata::onpe0() )
394 {
395 cout << " Wavefunction::set_deltaspin: nel+nempty too small" << endl;
396 cout << " Wavefunction::set_deltaspin: cannot set deltaspin" << endl;
397 return;
398 }
399 }
400 deallocate();
401 deltaspin_ = deltaspin;
402 compute_nst();
403 allocate();
404 resize();
405 init();
406 update_occ(0.0);
407 }
408
409 ////////////////////////////////////////////////////////////////////////////////
add_kpoint(D3vector kpoint,double weight)410 void Wavefunction::add_kpoint(D3vector kpoint, double weight)
411 {
412 for ( int i = 0; i < kpoint_.size(); i++ )
413 {
414 if ( length(kpoint - kpoint_[i]) < 1.e-6 )
415 {
416 if ( MPIdata::onpe0() )
417 cout << " Wavefunction::add_kpoint: kpoint already defined"
418 << endl;
419 return;
420 }
421 }
422
423 kpoint_.push_back(kpoint);
424 weight_.push_back(weight);
425
426 // Add SlaterDet
427 // determine on which block the SlaterDet should be added
428 const int ikp_new = kpoint_.size() - 1;
429 sd_.resize(nsp_loc());
430 const int ikp_loc = ikp_local(ikp_new);
431 if ( ikp_loc >= 0 )
432 {
433 // new kpoint is local to the current block
434 nkp_loc_[MPIdata::ikpb()]++;
435 ikp_global_.push_back(ikp_new);
436 for ( int isp_loc = 0; isp_loc < nsp_loc(); ++isp_loc )
437 {
438 sd_[isp_loc].push_back(new SlaterDet(sd_ctxt_,kpoint_[ikp_new]));
439 sd_[isp_loc][ikp_loc]->resize(cell_,refcell_,ecut_,
440 nst_[isp_global(isp_loc)]);
441 }
442
443 if ( nspin_ == 1 )
444 {
445 if ( nsp_loc() > 0 )
446 sd_[0][ikp_loc]->update_occ(nel_,nspin_);
447 }
448 else if ( nspin_ == 2 )
449 {
450 const int nocc_up = (nel_+1)/2+deltaspin_;
451 const int nocc_dn = nel_/2 - deltaspin_;
452 if ( nsp_loc() > 0 )
453 sd_[0][ikp_loc]->update_occ(nocc_up,nspin_);
454 if ( nsp_loc() > 1 )
455 sd_[1][ikp_loc]->update_occ(nocc_dn,nspin_);
456 }
457 else
458 {
459 // incorrect value of nspin_
460 assert(false);
461 }
462 }
463 }
464
465 ////////////////////////////////////////////////////////////////////////////////
del_kpoint(D3vector kpoint)466 void Wavefunction::del_kpoint(D3vector kpoint)
467 {
468 bool found = false;
469 vector<D3vector>::iterator pk = kpoint_.begin();
470 vector<double>::iterator pw = weight_.begin();
471 while ( !found && pk != kpoint_.end() )
472 {
473 if ( length(kpoint - *pk) < 1.e-6 )
474 {
475 found = true;
476 }
477 else
478 {
479 pk++;
480 pw++;
481 }
482 }
483 if ( !found )
484 {
485 if ( MPIdata::onpe0() )
486 cout << " Wavefunction::del_kpoint: no such kpoint"
487 << endl;
488 return;
489 }
490 deallocate();
491 kpoint_.erase(pk);
492 weight_.erase(pw);
493 allocate();
494 resize();
495 init();
496 update_occ(0.0);
497 }
498
499 ////////////////////////////////////////////////////////////////////////////////
move_kpoint(D3vector kpoint,D3vector newkpoint)500 void Wavefunction::move_kpoint(D3vector kpoint, D3vector newkpoint)
501 {
502 bool found = false;
503 int ikp = 0;
504 while ( !found && ikp < kpoint_.size() )
505 {
506 if ( length(kpoint_[ikp] - kpoint) < 1.e-6 )
507 {
508 found = true;
509 }
510 else
511 {
512 ikp++;
513 }
514 }
515 if ( !found )
516 {
517 if ( MPIdata::onpe0() )
518 cout << " Wavefunction::move_kpoint: no such kpoint"
519 << endl;
520 return;
521 }
522 // check if new kpoint position coincides with existing kpoint
523 for ( int i = 0; i < kpoint_.size(); i++ )
524 {
525 if ( length(newkpoint - kpoint_[i]) < 1.e-6 )
526 {
527 if ( MPIdata::onpe0() )
528 cout << " Wavefunction::move_kpoint: kpoint already defined "
529 << "at kpoint new position"
530 << endl;
531 return;
532 }
533 }
534
535 // ikp: global index of kpoint to be moved
536 const int ikp_loc = ikp_local(ikp);
537
538 // copy wavefunctions from old SlaterDet at kpoint to new SlaterDet
539 // at newkpoint
540 for ( int ispin = 0; ispin < nspin_; ispin++ )
541 {
542 const int isp_loc = isp_local(ispin);
543 if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
544 {
545 // this task holds kpoint ikp
546
547 // create new SlaterDet at newkpoint
548 SlaterDet *sd = sd_[isp_loc][ikp_loc];
549 SlaterDet *sdn = new SlaterDet(sd->context(),newkpoint);
550 sdn->resize(cell_,refcell_,ecut_,nst_[ispin]);
551 sdn->init();
552 // copy wave functions from old to new SlaterDet
553 const Basis& basis = sd_[isp_loc][ikp_loc]->basis();
554 const Basis& newbasis = sdn->basis();
555 // transform all states to real space and interpolate
556 int np0 = max(basis.np(0),newbasis.np(0));
557 int np1 = max(basis.np(1),newbasis.np(1));
558 int np2 = max(basis.np(2),newbasis.np(2));
559 FourierTransform ft1(basis,np0,np1,np2);
560 FourierTransform ft2(newbasis,np0,np1,np2);
561 // allocate real-space grid
562 valarray<complex<double> > tmpr(ft1.np012loc());
563 for ( int n = 0; n < sd->nstloc(); n++ )
564 {
565 for ( int i = 0; i < tmpr.size(); i++ )
566 tmpr[i] = 0.0;
567 ComplexMatrix& c = sd->c();
568 ComplexMatrix& cn = sdn->c();
569 ft1.backward(c.valptr(n*c.mloc()),&tmpr[0]);
570 // if the new kpoint is Gamma, remove the phase of the wavefunction
571 // using the phase of the first element of the array
572 if ( newbasis.real() )
573 {
574 double arg0 = arg(tmpr[0]);
575 // broadcast argument found on task 0
576 MPI_Bcast(&arg0,1,MPI_DOUBLE,0,basis.comm());
577 complex<double> z = polar(1.0,-arg0);
578 for ( int i = 0; i < tmpr.size(); i++ )
579 tmpr[i] *= z;
580 }
581 ft2.forward(&tmpr[0], cn.valptr(n*cn.mloc()));
582 }
583
584 // delete old SlaterDet
585 delete sd_[isp_loc][ikp_loc];
586 // reassign pointer
587 sd_[isp_loc][ikp_loc] = sdn;
588
589 // randomize new SlaterDet to break symmetry
590 sd_[isp_loc][ikp_loc]->randomize(1.e-4);
591
592 } // if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
593 } // for ispin
594
595 kpoint_[ikp] = newkpoint;
596
597 // update occupation numbers
598 if ( nspin_ == 1 )
599 {
600 if ( ( isp_local(0) >= 0 ) && ( ikp_loc >= 0 ) )
601 {
602 sd_[isp_local(0)][ikp_loc]->update_occ(nel_,nspin_);
603 }
604 }
605 else if ( nspin_ == 2 )
606 {
607 const int nocc_up = (nel_+1)/2+deltaspin_;
608 const int nocc_dn = nel_/2 - deltaspin_;
609
610 if ( ( isp_local(0) >= 0 ) && ( ikp_loc >= 0 ) )
611 sd_[isp_local(0)][ikp_loc]->update_occ(nocc_up,nspin_);
612
613 if ( ( isp_local(1) >= 0 ) && ( ikp_loc >= 0 ) )
614 sd_[isp_local(1)][ikp_loc]->update_occ(nocc_dn,nspin_);
615 }
616 else
617 {
618 // incorrect value of nspin_
619 assert(false);
620 }
621 }
622
623 ////////////////////////////////////////////////////////////////////////////////
randomize(double amplitude)624 void Wavefunction::randomize(double amplitude)
625 {
626 srand48((long int) MPIdata::rank());
627 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
628 {
629 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
630 {
631 sd_[isp_loc][ikp_loc]->randomize(amplitude);
632 }
633 }
634 }
635
636 ////////////////////////////////////////////////////////////////////////////////
update_occ(double temp)637 void Wavefunction::update_occ(double temp)
638 {
639 // update occupation numbers using eigenvalues in SlaterDet
640 if ( temp == 0.0 )
641 {
642 // zero temperature
643 if ( nspin_ == 1 )
644 {
645 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
646 {
647 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
648 {
649 sd_[isp_loc][ikp_loc]->update_occ(nel_,nspin_);
650 }
651 }
652 }
653 else if ( nspin_ == 2 )
654 {
655 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
656 {
657 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
658 {
659 if ( isp_global_[isp_loc] == 0 )
660 {
661 const int nocc_up = (nel_+1)/2+deltaspin_;
662 sd_[isp_loc][ikp_loc]->update_occ(nocc_up,nspin_);
663 }
664 else
665 {
666 const int nocc_dn = nel_/2 - deltaspin_;
667 sd_[isp_loc][ikp_loc]->update_occ(nocc_dn,nspin_);
668 }
669 }
670 }
671 }
672 else
673 {
674 // incorrect value of nspin_
675 assert(false);
676 }
677 }
678 else
679 {
680 // finite temperature
681 const double eVolt = 0.036749023; // 1 eV in Hartree
682 const int maxiter = 500;
683
684 // loop to find value of mu
685 double mu = 0.0;
686 double dmu = 2.0 * eVolt;
687 const double totalcharge = (double) nel_;
688 enum direction { up, down };
689 direction dir = up;
690
691 double rhosum = 0.0;
692 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
693 {
694 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
695 {
696 sd_[isp_loc][ikp_loc]->update_occ(nspin_,mu,temp);
697 const int ikpg = ikp_global_[ikp_loc];
698 rhosum += weight_[ikpg] * sd_[isp_loc][ikp_loc]->total_charge();
699 }
700 }
701 // sum over spin and kpoints
702 double tmpsum = 0.0;
703 MPI_Allreduce(&rhosum,&tmpsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::kp_sp_comm());
704 rhosum = tmpsum;
705
706 int niter = 0;
707 while ( niter < maxiter && fabs(rhosum - nel_) > 1.e-10 )
708 {
709 niter++;
710 if ( rhosum < totalcharge )
711 {
712 if ( dir == down ) dmu /= 2.0;
713 mu += dmu;
714 dir = up;
715 }
716 else
717 {
718 if ( dir == up ) dmu /= 2.0;
719 mu -= dmu;
720 dir = down;
721 }
722 rhosum = 0.0;
723 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
724 {
725 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
726 {
727 sd_[isp_loc][ikp_loc]->update_occ(nspin_,mu,temp);
728 int ikp = ikp_global_[ikp_loc];
729 rhosum += weight_[ikp] * sd_[isp_loc][ikp_loc]->total_charge();
730 }
731 }
732 MPI_Allreduce(&rhosum,&tmpsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::kp_sp_comm());
733 rhosum = tmpsum;
734 }
735
736 if ( niter == maxiter )
737 {
738 cout << "Wavefunction::update_occ: mu did not converge in "
739 << maxiter << " iterations" << endl;
740 MPI_Abort(MPIdata::comm(),1);
741 }
742
743 if ( MPIdata::onpe0() )
744 {
745 cout << " Wavefunction::update_occ: sum = "
746 << rhosum << endl;
747 cout << "<mu> " << setprecision(6) << mu / eVolt << " </mu>" << endl;
748 cout << "<occ_set>" << endl;
749 }
750
751 for ( int ispin = 0; ispin < nspin(); ++ispin )
752 {
753 const int isp_loc = isp_local(ispin);
754 for ( int ikp = 0; ikp < nkp(); ++ikp )
755 {
756 const int ikp_loc = ikp_local(ikp);
757 ostringstream ostr;
758 ostr.setf(ios::fixed,ios::floatfield);
759 ostr.setf(ios::right,ios::adjustfield);
760 int isrc = -1;
761 if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
762 {
763 if ( MPIdata::sd_rank() == 0 )
764 {
765 ostr.str("");
766 isrc = MPIdata::rank();
767 const int nst = sd(isp_loc,ikp_loc)->nst();
768 ostr << " <occ spin=\"" << ispin
769 << "\" kpoint=\""
770 << setprecision(8)
771 << sd(isp_loc,ikp_loc)->kpoint()
772 << "\" weight=\""
773 << setprecision(8)
774 << weight(ikp)
775 << "\" n=\"" << nst << "\">" << endl;
776 for ( int i = 0; i < nst; i++ )
777 {
778 ostr << setw(7) << setprecision(4)
779 << sd(isp_loc,ikp_loc)->occ(i);
780 if ( i%10 == 9 ) ostr << endl;
781 }
782 if ( nst%10 != 0 ) ostr << endl;
783 ostr << " </occ>" << endl;
784 }
785 }
786 cout0(ostr.str(),isrc);
787 MPI_Barrier(MPIdata::comm());
788 }
789 }
790 if ( MPIdata::onpe0() )
791 cout << "</occ_set>" << endl;
792 MPI_Barrier(MPIdata::comm());
793 }
794 }
795
796 ////////////////////////////////////////////////////////////////////////////////
gram(void)797 void Wavefunction::gram(void)
798 {
799 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
800 {
801 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
802 {
803 sd_[isp_loc][ikp_loc]->gram();
804 }
805 }
806 }
807
808 ////////////////////////////////////////////////////////////////////////////////
riccati(Wavefunction & wf)809 void Wavefunction::riccati(Wavefunction& wf)
810 {
811 assert(wf.sd_context() == sd_ctxt_);
812 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
813 {
814 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
815 {
816 sd_[isp_loc][ikp_loc]->riccati(*wf.sd_[isp_loc][ikp_loc]);
817 }
818 }
819 }
820
821 ////////////////////////////////////////////////////////////////////////////////
align(Wavefunction & wf)822 void Wavefunction::align(Wavefunction& wf)
823 {
824 assert(wf.sd_context() == sd_ctxt_);
825 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
826 {
827 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
828 {
829 sd_[isp_loc][ikp_loc]->align(*wf.sd_[isp_loc][ikp_loc]);
830 }
831 }
832 }
833
834 ////////////////////////////////////////////////////////////////////////////////
dot(const Wavefunction & wf) const835 complex<double> Wavefunction::dot(const Wavefunction& wf) const
836 {
837 assert(wf.sd_context() == sd_ctxt_);
838 complex<double> sum = 0.0;
839 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
840 {
841 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
842 {
843 const int ikpg = ikp_global(ikp_loc);
844 sum += weight_[ikpg] *
845 sd_[isp_loc][ikp_loc]->dot(*wf.sd_[isp_loc][ikp_loc]);
846 }
847 }
848 // sum over kpoint and spin comms
849 complex<double> tsum;
850 MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE_COMPLEX,MPI_SUM,MPIdata::kp_sp_comm());
851 return tsum;
852 }
853
854 ////////////////////////////////////////////////////////////////////////////////
diag(Wavefunction & dwf,bool eigvec)855 void Wavefunction::diag(Wavefunction& dwf, bool eigvec)
856 {
857 // subspace diagonalization of <*this | dwf>
858 // if eigvec==true, eigenvectors are computed and stored in *this, dwf is
859 // overwritten
860 for ( int isp_loc = 0; isp_loc < nsp_loc(); ++isp_loc )
861 {
862 assert( nst_[isp_global(isp_loc)] > 0 );
863 for ( int ikp_loc = 0; ikp_loc < nkp_loc(); ++ikp_loc )
864 {
865 // compute eigenvalues
866 if ( sd(isp_loc,ikp_loc)->basis().real() )
867 {
868 // proxy real matrices c, cp
869 DoubleMatrix c(sd(isp_loc,ikp_loc)->c());
870 DoubleMatrix cp(dwf.sd(isp_loc,ikp_loc)->c());
871
872 DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
873
874 // factor 2.0 in next line: G and -G
875 h.gemm('t','n',2.0,c,cp,0.0);
876 // rank-1 update correction
877 h.ger(-1.0,c,0,cp,0);
878
879 // cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
880 // << endl;
881 // cout << h;
882
883 #if 1
884 valarray<double> w(h.m());
885 if ( eigvec )
886 {
887 DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
888 h.syev('l',w,z);
889 //h.syevx('l',w,z,1.e-6);
890 cp = c;
891 c.gemm('n','n',1.0,cp,z,0.0);
892 }
893 else
894 {
895 h.syev('l',w);
896 }
897 #else
898 vector<double> w(h.m());
899 DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
900 const int maxsweep = 30;
901 int nsweep = jacobi(maxsweep,1.e-6,h,z,w);
902 if ( eigvec )
903 {
904 cp = c;
905 c.gemm('n','n',1.0,cp,z,0.0);
906 }
907 #endif
908 // set eigenvalues in SlaterDet
909 sd(isp_loc,ikp_loc)->set_eig(w);
910 }
911 else
912 {
913 ComplexMatrix& c(sd(isp_loc,ikp_loc)->c());
914 ComplexMatrix& cp(dwf.sd(isp_loc,ikp_loc)->c());
915 ComplexMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
916 h.gemm('c','n',1.0,c,cp,0.0);
917 // cout << " Hamiltonian at k = "
918 // << sd(ispin,ikp)->kpoint() << endl;
919 // cout << h;
920 valarray<double> w(h.m());
921 if ( eigvec )
922 {
923 ComplexMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
924 h.heev('l',w,z);
925 cp = c;
926 c.gemm('n','n',1.0,cp,z,0.0);
927 }
928 else
929 {
930 h.heev('l',w);
931 }
932 // set eigenvalues in SlaterDet
933 sd(isp_loc,ikp_loc)->set_eig(w);
934 }
935 }
936 }
937 }
938
939 ////////////////////////////////////////////////////////////////////////////////
print(ostream & os,string encoding,string tag) const940 void Wavefunction::print(ostream& os, string encoding, string tag) const
941 {
942 if ( MPIdata::onpe0() )
943 {
944 os << "<" << tag << " ecut=\"" << ecut_ << "\""
945 << " nspin=\"" << nspin_ << "\""
946 << " nel=\"" << nel_ << "\""
947 << " nempty=\"" << nempty_ << "\">" << endl;
948 os << setprecision(10);
949 os << "<domain a=\""
950 << cell_.a(0) << "\"\n b=\""
951 << cell_.a(1) << "\"\n c=\""
952 << cell_.a(2) << "\"/>" << endl;
953 os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
954 << " ny=\"" << sd_[0][0]->basis().np(1) << "\""
955 << " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
956 }
957
958 for ( int ispin = 0; ispin < nspin(); ++ispin )
959 {
960 const int isp_loc = isp_local(ispin);
961 for ( int ikp = 0; ikp < nkp(); ++ikp )
962 {
963 const int ikp_loc = ikp_local(ikp);
964 MPI_Barrier(MPIdata::comm());
965 if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
966 {
967 sd_[isp_loc][ikp_loc]->print(os,encoding,weight_[ikp],ispin,nspin_);
968 }
969 }
970 }
971 MPI_Barrier(MPIdata::comm());
972
973 if ( MPIdata::onpe0() )
974 os << "</" << tag << ">" << endl;
975 }
976
977 ///////////////////////////////////////////////////////////////////////////////
write(SharedFilePtr & sfp,string encoding,string tag) const978 void Wavefunction::write(SharedFilePtr& sfp, string encoding, string tag) const
979 {
980 assert(sizeof(size_t)==sizeof(MPI_Offset));
981 assert(sizeof(long long int)==sizeof(MPI_Offset));
982
983 sfp.sync();
984
985 if ( MPIdata::onpe0() )
986 {
987 ostringstream os;
988 os << "<" << tag << " ecut=\"" << ecut_ << "\""
989 << " nspin=\"" << nspin_ << "\""
990 << " nel=\"" << nel_ << "\""
991 << " nempty=\"" << nempty_ << "\">" << endl;
992 os << setprecision(10);
993 os << "<domain a=\""
994 << cell_.a(0) << "\"\n b=\""
995 << cell_.a(1) << "\"\n c=\""
996 << cell_.a(2) << "\"/>" << endl;
997 if ( refcell_.volume() != 0.0 )
998 {
999 os << "<reference_domain a=\""
1000 << refcell_.a(0) << "\"\n b=\""
1001 << refcell_.a(1) << "\"\n c=\""
1002 << refcell_.a(2) << "\"/>" << endl;
1003 }
1004 os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
1005 << " ny=\"" << sd_[0][0]->basis().np(1) << "\""
1006 << " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
1007 string str(os.str());
1008 size_t len = str.size();
1009 MPI_Status status;
1010 int err = MPI_File_write_at(sfp.file(),sfp.offset(),(void*)str.c_str(),
1011 len,MPI_CHAR,&status);
1012 if ( err != 0 )
1013 cout << " Wavefunction::write: error in MPI_File_write" << endl;
1014 sfp.advance(len);
1015 }
1016
1017 sfp.sync();
1018
1019 vector<vector<string> > sdstr;
1020 sdstr.resize(nsp_loc());
1021 for ( int isp_loc = 0; isp_loc < nsp_loc(); ++isp_loc )
1022 {
1023 const int ispin = isp_global(isp_loc);
1024 sdstr[isp_loc].resize(nkp_loc());
1025 for ( int ikp_loc = 0; ikp_loc < nkp_loc(); ++ikp_loc )
1026 {
1027 const int ikp = ikp_global(ikp_loc);
1028 // serialize sd[isp_loc][ikp_loc] into sdstr[isp_loc][ikp_loc]
1029 sd_[isp_loc][ikp_loc]->str(sdstr[isp_loc][ikp_loc],encoding,
1030 weight_[ikp],ispin,nspin_);
1031 }
1032 }
1033
1034 // sdstr now contains all data to be written
1035 // write data in order of increasing ispin, ikp
1036 for ( int ispin = 0; ispin < nspin(); ++ispin )
1037 {
1038 for ( int ikp = 0; ikp < nkp(); ++ikp )
1039 {
1040 MPI_Barrier(MPIdata::comm());
1041 const int isp_loc = isp_local(ispin);
1042 const int ikp_loc = ikp_local(ikp);
1043 const char* wbuf = 0;
1044 size_t len = 0;
1045 MPI_Offset off = 0;
1046 if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
1047 {
1048 wbuf = sdstr[isp_loc][ikp_loc].c_str();
1049 len = sdstr[isp_loc][ikp_loc].size();
1050 // compute offset of data on current task
1051 long long int local_offset = 0;
1052 long long int local_size = len;
1053 MPI_Scan(&local_size, &local_offset, 1,
1054 MPI_LONG_LONG, MPI_SUM, MPIdata::sd_comm());
1055 // correct for inclusive scan by subtracting local_size
1056 local_offset -= local_size;
1057 off = sfp.offset() + local_offset;
1058 }
1059 MPI_Status status;
1060 int err = MPI_File_write_at_all(sfp.file(),off,wbuf,len,
1061 MPI_CHAR,&status);
1062 int count;
1063 MPI_Get_count(&status,MPI_CHAR,&count);
1064 assert(count==len);
1065
1066 if ( err != 0 )
1067 cout << " Wavefunction::write: error in MPI_File_write_at_all" << endl;
1068 sfp.set_offset(off+len);
1069 sfp.sync();
1070 }
1071 }
1072
1073 if ( MPIdata::onpe0() )
1074 {
1075 ostringstream os;
1076 os << "</" << tag << ">" << endl;
1077 string str(os.str());
1078 size_t len = str.size();
1079 MPI_Status status;
1080 int err = MPI_File_write_at(sfp.file(),sfp.offset(),(void*)str.c_str(),
1081 len,MPI_CHAR,&status);
1082 if ( err != 0 )
1083 cout << " Wavefunction::write: error in MPI_File_write" << endl;
1084 sfp.advance(len);
1085 }
1086 }
1087
1088 ////////////////////////////////////////////////////////////////////////////////
info(ostream & os,string tag) const1089 void Wavefunction::info(ostream& os, string tag) const
1090 {
1091 if ( MPIdata::onpe0() )
1092 {
1093 os << "<" << tag << " ecut=\"" << ecut_ << "\""
1094 << " nspin=\"" << nspin_ << "\""
1095 << " nel=\"" << nel_ << "\""
1096 << " nempty=\"" << nempty_ << "\">" << endl;
1097 os.setf(ios::fixed,ios::floatfield);
1098 os << "<cell a=\""
1099 << setprecision(6) << cell_.a(0) << "\"\n b=\""
1100 << cell_.a(1) << "\"\n c=\""
1101 << cell_.a(2) << "\"/>" << endl;
1102 os << " reciprocal lattice vectors" << endl
1103 << setprecision(6)
1104 << " " << cell_.b(0) << endl
1105 << " " << cell_.b(1) << endl
1106 << " " << cell_.b(2) << endl;
1107 os << "<refcell a=\""
1108 << refcell_.a(0) << "\"\n b=\""
1109 << refcell_.a(1) << "\"\n c=\""
1110 << refcell_.a(2) << "\"/>" << endl;
1111 os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
1112 << " ny=\"" << sd_[0][0]->basis().np(1) << "\""
1113 << " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
1114 }
1115
1116 for ( int ispin = 0; ispin < nspin(); ++ispin )
1117 {
1118 const int isp_loc = isp_local(ispin);
1119 for ( int ikp = 0; ikp < nkp(); ++ikp )
1120 {
1121 const int ikp_loc = ikp_local(ikp);
1122 string s;
1123 int isrc = -1;
1124 if ( ( isp_loc >= 0 ) && ( ikp_loc >= 0 ) )
1125 {
1126 s = sd_[isp_loc][ikp_loc]->info();
1127 if ( MPIdata::sd_rank() == 0 )
1128 isrc = MPIdata::rank();
1129 }
1130 cout0(s,isrc);
1131 MPI_Barrier(MPIdata::comm());
1132 }
1133 }
1134
1135 if ( MPIdata::onpe0() )
1136 os << "</" << tag << ">" << endl;
1137 }
1138
1139 ////////////////////////////////////////////////////////////////////////////////
operator <<(ostream & os,const Wavefunction & wf)1140 ostream& operator<<(ostream& os, const Wavefunction& wf)
1141 {
1142 wf.print(os,"text","wavefunction");
1143 return os;
1144 }
1145
1146 ////////////////////////////////////////////////////////////////////////////////
operator =(const Wavefunction & wf)1147 Wavefunction& Wavefunction::operator=(const Wavefunction& wf)
1148 {
1149 if ( this == &wf ) return *this;
1150 assert(sd_ctxt_ == wf.sd_ctxt_);
1151 assert(nel_ == wf.nel_);
1152 assert(nempty_== wf.nempty_);
1153 assert(nspin_ == wf.nspin_);
1154 assert(deltaspin_ == wf.deltaspin_);
1155 assert(refcell_ == wf.refcell_);
1156 assert(ecut_ == wf.ecut_);
1157
1158 cell_ = wf.cell_;
1159 weight_ = wf.weight_;
1160 kpoint_ = wf.kpoint_;
1161
1162 for ( int isp_loc = 0; isp_loc < sd_.size(); ++isp_loc )
1163 {
1164 for ( int ikp_loc = 0; ikp_loc < sd_[isp_loc].size(); ++ikp_loc )
1165 {
1166 *sd_[isp_loc][ikp_loc] = *wf.sd_[isp_loc][ikp_loc];
1167 }
1168 }
1169 return *this;
1170 }
1171