Wavefunction.C 21.4 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// Wavefunction.C
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: Wavefunction.C,v 1.32 2008-04-15 01:36:44 fgygi Exp $
Francois Gygi committed
7 8 9

#include "Wavefunction.h"
#include "SlaterDet.h"
10
#include "jacobi.h"
11
#include "SharedFilePtr.h"
Francois Gygi committed
12 13 14 15 16 17
#include <vector>
#include <iomanip>
#include <sstream>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
18 19
Wavefunction::Wavefunction(const Context& ctxt) : ctxt_(ctxt), nel_(0),
nempty_(0), nspin_(1), deltaspin_(0), ecut_(0.0), nrowmax_(32)
Francois Gygi committed
20 21 22 23 24 25 26
{
  // create a default wavefunction: one k point, k=0
  kpoint_.resize(1);
  kpoint_[0] = D3vector(0,0,0);
  weight_.resize(1);
  weight_[0] = 1.0;
  compute_nst();
27
  create_contexts();
Francois Gygi committed
28 29 30 31
  allocate();
}

////////////////////////////////////////////////////////////////////////////////
32 33 34
Wavefunction::Wavefunction(const Wavefunction& wf) : ctxt_(wf.ctxt_),
nel_(wf.nel_), nempty_(wf.nempty_), nspin_(wf.nspin_),
deltaspin_(wf.deltaspin_), nrowmax_(wf.nrowmax_),
35
cell_(wf.cell_), refcell_(wf.refcell_),
Francois Gygi committed
36 37
ecut_(wf.ecut_), weight_(wf.weight_), kpoint_(wf.kpoint_)
{
38
  // Create a Wavefunction using the dimensions of the argument
39

Francois Gygi committed
40
  compute_nst();
41 42

  // Next lines: do special allocation of contexts to ensure that
Francois Gygi committed
43
  // contexts are same as those of wf
44
  spincontext_ = new Context(*wf.spincontext_);
45 46
  kpcontext_ = new Context(*wf.kpcontext_);
  sdcontext_ = new Context(*wf.sdcontext_);
Francois Gygi committed
47

48
  allocate();
49

50 51
  resize(cell_,refcell_,ecut_);
  reset();
Francois Gygi committed
52 53 54 55 56 57
}

////////////////////////////////////////////////////////////////////////////////
Wavefunction::~Wavefunction()
{
  deallocate();
58 59 60
  delete spincontext_;
  delete kpcontext_;
  delete sdcontext_;
Francois Gygi committed
61 62 63 64 65
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::allocate(void)
{
66
  // create SlaterDets using kpoint list
Francois Gygi committed
67 68
  const int nkp = kpoint_.size();
  assert(nspin_ == 1);
Francois Gygi committed
69 70 71 72 73
  sd_.resize(nspin_);
  sd_[0].resize(nkp);

  for ( int ikp = 0; ikp < nkp; ikp++ )
    sd_[0][ikp] = new SlaterDet(*sdcontext_,kpoint_[ikp]);
Francois Gygi committed
74 75 76 77 78
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::deallocate(void)
{
79 80
  assert (nspin_==1);
  for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
81
  {
82
    delete sd_[0][ikp];
Francois Gygi committed
83
  }
84
  sd_[0].resize(0);
Francois Gygi committed
85 86 87 88 89 90 91 92 93
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::clear(void)
{
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < nkp(); ikp++ )
    {
Francois Gygi committed
94
      sd(ispin,ikp)->c().clear();
Francois Gygi committed
95 96 97
    }
  }
}
Francois Gygi committed
98

Francois Gygi committed
99
////////////////////////////////////////////////////////////////////////////////
100
int Wavefunction::nkp(void) const { return kpoint_.size(); }
Francois Gygi committed
101 102 103 104 105 106

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nel() const { return nel_; } // total number of electrons

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nst() const
107
{
Francois Gygi committed
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
  if ( nspin_ == 1 )
    return nst_[0];
  else
    return nst_[0]+nst_[1];
}

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nst(int ispin) const
{
  assert(ispin >= 0 && ispin < 2);
  return nst_[ispin];
}

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nempty() const { return nempty_; } // number of empty states

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
125
int Wavefunction::nspin() const { return nspin_; } // number of spin components
Francois Gygi committed
126 127

////////////////////////////////////////////////////////////////////////////////
128
double Wavefunction::entropy(void) const
Francois Gygi committed
129
{
Francois Gygi committed
130 131 132 133 134
  double sum = 0.0;
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < nkp(); ikp++ )
    {
135
      sum += weight_[ikp] * sd(ispin,ikp)->entropy(nspin_);
Francois Gygi committed
136 137 138
    }
  }
  return sum;
Francois Gygi committed
139
}
Francois Gygi committed
140 141

////////////////////////////////////////////////////////////////////////////////
142
void Wavefunction::resize(const UnitCell& cell, const UnitCell& refcell,
Francois Gygi committed
143 144
  double ecut)
{
145 146 147 148 149 150 151
  try
  {
    // resize all SlaterDets using cell, refcell, ecut and nst_[ispin]
    for ( int ispin = 0; ispin < nspin_; ispin++ )
    {
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
      {
Francois Gygi committed
152
        sd_[ispin][ikp]->resize(cell,refcell,ecut,nst_[ispin]);
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
      }
    }
    cell_ = cell;
    refcell_ = refcell;
    ecut_ = ecut;
  }
  catch ( const SlaterDetException& sdex )
  {
    cout << " Wavefunction: SlaterDetException during resize: " << endl
         << sdex.msg << endl;
    // no resize took place
    return;
  }
  catch ( bad_alloc )
  {
    cout << " Wavefunction: insufficient memory for resize operation" << endl;
    return;
  }
171

172
}
173

174 175 176 177 178 179 180 181
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::reset(void)
{
  // reset all SlaterDets
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
182
      sd_[ispin][ikp]->reset();
183 184
    }
  }
Francois Gygi committed
185
}
186

Francois Gygi committed
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::compute_nst(void)
{
  // recompute nst from nel_, deltaspin_, nempty_

  nst_.resize(nspin_);
  if ( nspin_ == 1 )
  {
    nst_[0] = ( nel_ + 1 ) / 2 + nempty_;
  }
  else
  {
    // nspin == 2
    nst_[0] = ( nel_ + 1 ) / 2 + deltaspin_ + nempty_;
    nst_[1] = nel_ / 2 - deltaspin_ + nempty_;
  }
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nel(int nel)
{
208
  if ( nel == nel_ ) return;
Francois Gygi committed
209 210 211 212 213
  if ( nel < 0 )
  {
    cout << " Wavefunction::set_nel: nel < 0" << endl;
    return;
  }
214

Francois Gygi committed
215 216
  nel_ = nel;
  compute_nst();
217
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
218
}
219

Francois Gygi committed
220 221 222
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nempty(int nempty)
{
223
  if ( nempty == nempty_ ) return;
Francois Gygi committed
224 225 226 227 228 229 230
  if ( nempty < 0 )
  {
    cout << " Wavefunction::set_nempty: negative value" << endl;
    return;
  }
  nempty_ = nempty;
  compute_nst();
231
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
232 233 234 235 236 237 238
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nspin(int nspin)
{
  assert(nspin==1 || nspin==2);
  if ( nspin == nspin_ ) return;
239

Francois Gygi committed
240 241
  deallocate();
  cout << " Wavefunction::set_nspin: " << nspin << " deallocate done" << endl;
242

Francois Gygi committed
243
  nspin_ = nspin;
244

245
  compute_nst();
Francois Gygi committed
246 247
  allocate();
  cout << " Wavefunction::set_nspin: " << nspin << " allocate done" << endl;
248 249
  resize(cell_,refcell_,ecut_);
  reset();
Francois Gygi committed
250 251
}

252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::create_contexts(void)
{
  // determine dimensions of sdcontext
  assert(nrowmax_>0);
  const int size = ctxt_.size();
  int npr = nrowmax_;
  while ( size%npr != 0 ) npr--;
  // npr now divides size
  int npc = size/npr;

  spincontext_ = new Context(npr,npc);
  kpcontext_ = new Context(*spincontext_);
  sdcontext_ = new Context(*kpcontext_);
}

Francois Gygi committed
268 269 270 271 272
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nrowmax(int n)
{
  if ( n > ctxt_.size() )
  {
273 274
    if ( ctxt_.onpe0() )
      cout << " Wavefunction::set_nrowmax: nrowmax > ctxt_.size()" << endl;
Francois Gygi committed
275 276
    return;
  }
277

Francois Gygi committed
278
  deallocate();
279 280 281
  delete spincontext_;
  delete kpcontext_;
  delete sdcontext_;
Francois Gygi committed
282
  nrowmax_ = n;
283
  create_contexts();
284
  compute_nst();
Francois Gygi committed
285
  allocate();
286 287
  resize(cell_,refcell_,ecut_);
  reset();
Francois Gygi committed
288
}
289

Francois Gygi committed
290 291 292 293 294
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::add_kpoint(D3vector kpoint, double weight)
{
  for ( int i = 0; i < kpoint_.size(); i++ )
  {
Francois Gygi committed
295
    if ( length(kpoint - kpoint_[i]) < 1.e-6 )
Francois Gygi committed
296
    {
Francois Gygi committed
297 298
      if ( ctxt_.onpe0() )
        cout << " Wavefunction::add_kpoint: kpoint already defined"
Francois Gygi committed
299
           << endl;
Francois Gygi committed
300
      return;
Francois Gygi committed
301 302
    }
  }
303

Francois Gygi committed
304 305 306
  deallocate();
  kpoint_.push_back(kpoint);
  weight_.push_back(weight);
307

Francois Gygi committed
308
  allocate();
309 310
  resize(cell_,refcell_,ecut_);
  reset();
Francois Gygi committed
311
  update_occ(0.0);
Francois Gygi committed
312 313 314 315 316
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::del_kpoint(D3vector kpoint)
{
Francois Gygi committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
  bool found = false;
  vector<D3vector>::iterator pk = kpoint_.begin();
  vector<double>::iterator pw = weight_.begin();
  while ( !found && pk != kpoint_.end() )
  {
    if ( length(kpoint - *pk) < 1.e-6 )
    {
      found = true;
    }
    else
    {
      pk++;
      pw++;
    }
  }
  if ( !found )
  {
    if ( ctxt_.onpe0() )
      cout << " Wavefunction::del_kpoint: no such kpoint"
         << endl;
    return;
  }
  deallocate();
  kpoint_.erase(pk);
  weight_.erase(pw);
  allocate();
  resize(cell_,refcell_,ecut_);
  reset();
  update_occ(0.0);
Francois Gygi committed
346 347 348 349 350 351 352 353 354
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::randomize(double amplitude)
{
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
355
      sd_[ispin][ikp]->randomize(amplitude);
Francois Gygi committed
356 357 358 359 360 361 362
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::update_occ(double temp)
{
Francois Gygi committed
363 364
  // update occupation numbers using eigenvalues in SlaterDet
  if ( temp == 0.0 )
Francois Gygi committed
365
  {
Francois Gygi committed
366 367
    // zero temperature
    if ( nspin_ == 1 )
Francois Gygi committed
368
    {
Francois Gygi committed
369
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
370
      {
Francois Gygi committed
371
        sd_[0][ikp]->update_occ(nel_,nspin_);
Francois Gygi committed
372 373
      }
    }
Francois Gygi committed
374
    else if ( nspin_ == 2 )
Francois Gygi committed
375
    {
Francois Gygi committed
376 377 378 379
      const int nocc_up = (nel_+1)/2+deltaspin_;
      const int nocc_dn = nel_/2 - deltaspin_;
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
      {
Francois Gygi committed
380 381
        sd_[0][ikp]->update_occ(nocc_up,nspin_);
        sd_[1][ikp]->update_occ(nocc_dn,nspin_);
Francois Gygi committed
382
      }
Francois Gygi committed
383 384 385
    }
    else
    {
Francois Gygi committed
386 387
      // incorrect value of nspin_
      assert(false);
Francois Gygi committed
388
    }
Francois Gygi committed
389 390 391 392 393
  }
  else
  {
    // finite temperature
    const double eVolt = 0.036749023; // 1 eV in Hartree
394
    const int maxiter = 500;
395

Francois Gygi committed
396 397
    // loop to find value of mu
    double mu = 0.0;
398
    double dmu = 2.0 * eVolt;
Francois Gygi committed
399 400 401 402 403
    const double totalcharge = (double) nel_;
    enum direction { up, down };
    direction dir = up;

    double rhosum = 0.0;
Francois Gygi committed
404 405 406 407
    for ( int ispin = 0; ispin < nspin_; ispin++ )
    {
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
      {
Francois Gygi committed
408
        sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
409
        rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
410 411
      }
    }
412

Francois Gygi committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
    int niter = 0;
    while ( niter < maxiter && fabs(rhosum - nel_) > 1.e-10 )
    {
      niter++;
      if ( rhosum < totalcharge )
      {
        if ( dir == down ) dmu /= 2.0;
        mu += dmu;
        dir = up;
      }
      else
      {
        if ( dir == up ) dmu /= 2.0;
        mu -= dmu;
        dir = down;
      }
      rhosum = 0.0;
      for ( int ispin = 0; ispin < nspin_; ispin++ )
      {
        for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
        {
Francois Gygi committed
434
          sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
435
          rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
436 437 438
        }
      }
    }
439

Francois Gygi committed
440 441 442 443 444 445
    if ( niter == maxiter )
    {
      cout << "Wavefunction::update_occ: mu did not converge in "
           << maxiter << " iterations" << endl;
      ctxt_.abort(1);
    }
Francois Gygi committed
446

Francois Gygi committed
447
    if ( ctxt_.onpe0() )
Francois Gygi committed
448
    {
449 450 451 452
      cout << " Wavefunction::update_occ: sum = "
           << rhosum << endl;
      cout << " Wavefunction::update_occ: mu = "
           << setprecision(4) << mu / eVolt << " eV" << endl;
Francois Gygi committed
453 454 455

      cout.setf(ios::right,ios::adjustfield);
      cout.setf(ios::fixed,ios::floatfield);
456

457
      cout << " Wavefunction::update_occ: occupation numbers" << endl;
Francois Gygi committed
458
      for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
459
      {
Francois Gygi committed
460
        for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
461
        {
462
          cout << " k = " << kpoint_[ikp] << endl;
Francois Gygi committed
463
          for ( int n = 0; n < sd_[ispin][ikp]->nst(); n++ )
Francois Gygi committed
464
          {
Francois Gygi committed
465 466
            cout << setw(7) << setprecision(4) << sd_[ispin][ikp]->occ(n);
            if ( ( n%10 ) == 9 ) cout << endl;
Francois Gygi committed
467
          }
468 469
          if ( sd_[ispin][ikp]->nst() % 10 != 0 )
            cout << endl;
Francois Gygi committed
470 471 472 473 474 475 476
        }
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
477
void Wavefunction::gram(void)
Francois Gygi committed
478
{
Francois Gygi committed
479
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
480 481 482
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
483
      sd_[ispin][ikp]->gram();
Francois Gygi committed
484 485
    }
  }
Francois Gygi committed
486 487 488 489 490 491 492
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::riccati(Wavefunction& wf)
{
  assert(wf.context() == ctxt_);
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
493 494 495
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
496
      sd_[ispin][ikp]->riccati(*wf.sd_[ispin][ikp]);
Francois Gygi committed
497 498 499 500 501
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
502
void Wavefunction::align(Wavefunction& wf)
Francois Gygi committed
503
{
Francois Gygi committed
504
  assert(wf.context() == ctxt_);
Francois Gygi committed
505 506 507 508
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
509
      sd_[ispin][ikp]->align(*wf.sd_[ispin][ikp]);
Francois Gygi committed
510 511 512 513 514
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
515
complex<double> Wavefunction::dot(const Wavefunction& wf) const
Francois Gygi committed
516 517
{
  assert(wf.context() == ctxt_);
Francois Gygi committed
518
  complex<double> sum = 0.0;
Francois Gygi committed
519 520 521 522
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
523
      sum += sd_[ispin][ikp]->dot(*wf.sd_[ispin][ikp]);
Francois Gygi committed
524 525 526 527 528 529 530 531 532
    }
  }
  return sum;
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::diag(Wavefunction& dwf, bool eigvec)
{
  // subspace diagonalization of <*this | dwf>
533
  // if eigvec==true, eigenvectors are computed and stored in *this, dwf is
Francois Gygi committed
534 535 536 537 538
  // overwritten
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < nkp(); ikp++ )
    {
Francois Gygi committed
539 540
      // compute eigenvalues
      if ( sd(ispin,ikp)->basis().real() )
Francois Gygi committed
541
      {
Francois Gygi committed
542 543 544
        // proxy real matrices c, cp
        DoubleMatrix c(sd(ispin,ikp)->c());
        DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
545

Francois Gygi committed
546
        DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
547

Francois Gygi committed
548 549 550 551
        // factor 2.0 in next line: G and -G
        h.gemm('t','n',2.0,c,cp,0.0);
        // rank-1 update correction
        h.ger(-1.0,c,0,cp,0);
552

Francois Gygi committed
553 554 555
        // cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
        //      << endl;
        // cout << h;
556

557
#if 1
Francois Gygi committed
558 559 560 561 562 563 564 565 566 567 568 569 570
        valarray<double> w(h.m());
        if ( eigvec )
        {
          DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
          h.syev('l',w,z);
          //h.syevx('l',w,z,1.e-6);
          cp = c;
          c.gemm('n','n',1.0,cp,z,0.0);
        }
        else
        {
          h.syev('l',w);
        }
571
#else
Francois Gygi committed
572 573 574 575 576 577 578 579 580
        vector<double> w(h.m());
        DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
        const int maxsweep = 30;
        int nsweep = jacobi(maxsweep,1.e-6,h,z,w);
        if ( eigvec )
        {
          cp = c;
          c.gemm('n','n',1.0,cp,z,0.0);
        }
581
#endif
Francois Gygi committed
582 583 584 585 586 587 588 589 590
        // set eigenvalues in SlaterDet
        sd(ispin,ikp)->set_eig(w);
      }
      else
      {
        ComplexMatrix& c(sd(ispin,ikp)->c());
        ComplexMatrix& cp(dwf.sd(ispin,ikp)->c());
        ComplexMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
        h.gemm('c','n',1.0,c,cp,0.0);
591 592 593
        // cout << " Hamiltonian at k = "
        //      << sd(ispin,ikp)->kpoint() << endl;
        // cout << h;
Francois Gygi committed
594
        valarray<double> w(h.m());
595
        if ( eigvec )
Francois Gygi committed
596
        {
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618
#if DEBUG
          ComplexMatrix hcopy(h);
#endif
          ComplexMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
          h.heev('l',w,z);
          cp = c;
          c.gemm('n','n',1.0,cp,z,0.0);
#if DEBUG
          // check that z contains eigenvectors of h
          // diagonal matrix with eigenvalues on diagonal
          ComplexMatrix d(c.context(),c.n(),c.n(),c.nb(),c.nb());
          // the following test works only on one task
          assert(ctxt_.size()==1);
          for ( int i = 0; i < d.m(); i++ )
            d[i+d.n()*i] = w[i];
          ComplexMatrix dz(c.context(),c.n(),c.n(),c.nb(),c.nb());
          dz.gemm('n','c',1.0,d,z,0.0);
          ComplexMatrix zdz(c.context(),c.n(),c.n(),c.nb(),c.nb());
          zdz.gemm('n','n',1.0,z,dz,0.0);
          // zdz should be equal to hcopy
          zdz -= hcopy;
          cout << " heev: norm of error: " << zdz.nrm2() << endl;
Francois Gygi committed
619
#endif
620 621 622 623 624
        }
        else
        {
          h.heev('l',w);
        }
Francois Gygi committed
625 626
        // set eigenvalues in SlaterDet
        sd(ispin,ikp)->set_eig(w);
Francois Gygi committed
627 628 629 630 631 632
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
633
void Wavefunction::print(ostream& os, string encoding, string tag) const
Francois Gygi committed
634
{
635 636 637 638 639 640
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
641
    os << setprecision(10);
642 643 644 645 646 647 648 649
    os << "<domain a=\""
       << cell_.a(0) << "\"\n        b=\""
       << cell_.a(1) << "\"\n        c=\""
       << cell_.a(2) << "\"/>" << endl;
    os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
       <<      " ny=\"" << sd_[0][0]->basis().np(1) << "\""
       <<      " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
  }
650

Francois Gygi committed
651 652 653
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
654
      sd_[ispin][ikp]->print(os,encoding,weight_[ikp],ispin,nspin_);
Francois Gygi committed
655
  }
656

657 658
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
659 660 661
}

////////////////////////////////////////////////////////////////////////////////
662
void Wavefunction::write(SharedFilePtr& sfp, string encoding, string tag) const
Francois Gygi committed
663
{
664
  sfp.sync();
665

Francois Gygi committed
666 667 668 669 670 671 672
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
673
    os << setprecision(10);
Francois Gygi committed
674 675 676 677
    os << "<domain a=\""
       << cell_.a(0) << "\"\n        b=\""
       << cell_.a(1) << "\"\n        c=\""
       << cell_.a(2) << "\"/>" << endl;
678 679 680 681 682 683 684
    if ( refcell_.volume() != 0.0 )
    {
      os << "<reference_domain a=\""
         << refcell_.a(0) << "\"\n        b=\""
         << refcell_.a(1) << "\"\n        c=\""
         << refcell_.a(2) << "\"/>" << endl;
    }
Francois Gygi committed
685 686 687 688
    os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
       <<      " ny=\"" << sd_[0][0]->basis().np(1) << "\""
       <<      " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
    string str(os.str());
689 690
    int len = str.size();
    MPI_Status status;
691 692
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
693
    if ( err != 0 )
694
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
695
    sfp.advance(len);
Francois Gygi committed
696 697 698 699 700
  }

  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
701
    {
702
      sd_[ispin][ikp]->write(sfp,encoding,weight_[ikp],ispin,nspin_);
703
    }
Francois Gygi committed
704
  }
705

706
  sfp.sync();
707

Francois Gygi committed
708 709 710 711 712
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "</" << tag << ">" << endl;
    string str(os.str());
713 714
    int len = str.size();
    MPI_Status status;
715 716
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
717
    if ( err != 0 )
718
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
719
    sfp.advance(len);
Francois Gygi committed
720 721 722 723
  }
}

////////////////////////////////////////////////////////////////////////////////
724
void Wavefunction::info(ostream& os, string tag) const
Francois Gygi committed
725
{
726 727 728 729 730 731
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
732
    os.setf(ios::fixed,ios::floatfield);
733
    os << "<cell a=\""
734 735
       << setprecision(6) << cell_.a(0) << "\"\n      b=\""
       << cell_.a(1) << "\"\n      c=\""
736
       << cell_.a(2) << "\"/>" << endl;
737
    os << " reciprocal lattice vectors" << endl
Francois Gygi committed
738 739 740
       << setprecision(6)
       << " " << cell_.b(0) << endl
       << " " << cell_.b(1) << endl
741
       << " " << cell_.b(2) << endl;
742
    os << "<refcell a=\""
743 744
       << refcell_.a(0) << "\"\n         b=\""
       << refcell_.a(1) << "\"\n         c=\""
745
       << refcell_.a(2) << "\"/>" << endl;
746 747 748 749
    os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
       <<      " ny=\"" << sd_[0][0]->basis().np(1) << "\""
       <<      " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
  }
750

Francois Gygi committed
751 752 753
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
754
    {
755
      if ( ctxt_.onpe0() )
756 757
        cout << " kpoint: " << kpoint_[ikp] << " weight: " << weight_[ikp]
             << endl;
Francois Gygi committed
758
      sd_[ispin][ikp]->info(os);
Francois Gygi committed
759
    }
Francois Gygi committed
760
  }
761

762 763
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
764 765 766
}

////////////////////////////////////////////////////////////////////////////////
767
ostream& operator<<(ostream& os, const Wavefunction& wf)
Francois Gygi committed
768 769 770 771 772 773 774 775 776 777 778 779 780
{
  wf.print(os,"text","wavefunction");
  return os;
}

////////////////////////////////////////////////////////////////////////////////
Wavefunction& Wavefunction::operator=(const Wavefunction& wf)
{
  if ( this == &wf ) return *this;
  assert(ctxt_ == wf.ctxt_);
  assert(nel_ == wf.nel_);
  assert(nempty_== wf.nempty_);
  assert(nspin_ == wf.nspin_);
781
  assert(nrowmax_ == wf.nrowmax_);
Francois Gygi committed
782 783 784 785
  assert(deltaspin_ == wf.deltaspin_);
  assert(cell_ == wf.cell_);
  assert(refcell_ == wf.refcell_);
  assert(ecut_ == wf.ecut_);
786

Francois Gygi committed
787 788
  weight_ = wf.weight_;
  kpoint_ = wf.kpoint_;
789

790 791 792 793
  assert(*spincontext_ == *wf.spincontext_);
  assert(*kpcontext_   == *wf.kpcontext_);
  assert(*sdcontext_   == *wf.sdcontext_);

Francois Gygi committed
794 795 796 797
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
798
      *sd_[ispin][ikp] = *wf.sd_[ispin][ikp];
Francois Gygi committed
799 800 801 802
    }
  }
  return *this;
}