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