Wavefunction.C 23.8 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20
// Wavefunction.C
//
////////////////////////////////////////////////////////////////////////////////

#include "Wavefunction.h"
#include "SlaterDet.h"
21
#include "jacobi.h"
22
#include "SharedFilePtr.h"
Francois Gygi committed
23 24 25 26 27 28
#include <vector>
#include <iomanip>
#include <sstream>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
29 30
Wavefunction::Wavefunction(const Context& ctxt) : ctxt_(ctxt), nel_(0),
nempty_(0), nspin_(1), deltaspin_(0), ecut_(0.0), nrowmax_(32)
Francois Gygi committed
31 32 33 34 35 36 37
{
  // 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();
38
  create_contexts();
Francois Gygi committed
39 40 41 42
  allocate();
}

////////////////////////////////////////////////////////////////////////////////
43 44 45
Wavefunction::Wavefunction(const Wavefunction& wf) : ctxt_(wf.ctxt_),
nel_(wf.nel_), nempty_(wf.nempty_), nspin_(wf.nspin_),
deltaspin_(wf.deltaspin_), nrowmax_(wf.nrowmax_),
46
cell_(wf.cell_), refcell_(wf.refcell_),
Francois Gygi committed
47 48
ecut_(wf.ecut_), weight_(wf.weight_), kpoint_(wf.kpoint_)
{
49
  // Create a Wavefunction using the dimensions of the argument
50

Francois Gygi committed
51
  compute_nst();
52 53

  // Next lines: do special allocation of contexts to ensure that
Francois Gygi committed
54
  // contexts are same as those of wf
55
  spincontext_ = new Context(*wf.spincontext_);
56 57
  kpcontext_ = new Context(*wf.kpcontext_);
  sdcontext_ = new Context(*wf.sdcontext_);
Francois Gygi committed
58

59
  allocate();
60

61
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
62
  init();
Francois Gygi committed
63 64 65 66 67 68
}

////////////////////////////////////////////////////////////////////////////////
Wavefunction::~Wavefunction()
{
  deallocate();
69 70 71
  delete spincontext_;
  delete kpcontext_;
  delete sdcontext_;
Francois Gygi committed
72 73 74 75 76
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::allocate(void)
{
77
  // create SlaterDets using kpoint list
Francois Gygi committed
78
  const int nkp = kpoint_.size();
Francois Gygi committed
79
  sd_.resize(nspin_);
80
  for ( int ispin = 0; ispin < nspin_; ispin++ )
81 82 83 84 85 86 87
  {
    sd_[ispin].resize(nkp);
    for ( int ikp = 0; ikp < nkp; ikp++ )
    {
      sd_[ispin][ikp] = new SlaterDet(*sdcontext_,kpoint_[ikp]);
    }
  }
Francois Gygi committed
88 89 90 91 92
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::deallocate(void)
{
93
  for ( int ispin = 0; ispin < sd_.size(); ispin++ )
Francois Gygi committed
94
  {
95 96 97 98 99
    for ( int ikp = 0; ikp < sd_[ispin].size(); ikp++ )
    {
      delete sd_[ispin][ikp];
    }
    sd_[ispin].resize(0);
Francois Gygi committed
100 101 102 103
  }
}

////////////////////////////////////////////////////////////////////////////////
104
int Wavefunction::nkp(void) const { return kpoint_.size(); }
Francois Gygi committed
105 106 107 108 109 110

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

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nst() const
111
{
Francois Gygi committed
112 113 114 115 116 117 118 119 120 121 122 123 124 125
  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];
}

////////////////////////////////////////////////////////////////////////////////
126 127 128 129
int Wavefunction::nempty() const { return nempty_; }

////////////////////////////////////////////////////////////////////////////////
int Wavefunction::nspin() const { return nspin_; }
Francois Gygi committed
130 131

////////////////////////////////////////////////////////////////////////////////
132
int Wavefunction::deltaspin() const { return deltaspin_; }
Francois Gygi committed
133 134

////////////////////////////////////////////////////////////////////////////////
135
double Wavefunction::entropy(void) const
Francois Gygi committed
136
{
Francois Gygi committed
137 138 139 140 141
  double sum = 0.0;
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < nkp(); ikp++ )
    {
142
      sum += weight_[ikp] * sd(ispin,ikp)->entropy(nspin_);
Francois Gygi committed
143 144 145
    }
  }
  return sum;
Francois Gygi committed
146
}
Francois Gygi committed
147 148

////////////////////////////////////////////////////////////////////////////////
149
void Wavefunction::resize(const UnitCell& cell, const UnitCell& refcell,
Francois Gygi committed
150 151
  double ecut)
{
152 153 154 155 156 157 158
  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
159
        sd_[ispin][ikp]->resize(cell,refcell,ecut,nst_[ispin]);
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
      }
    }
    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;
  }
178

179
}
180

181
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
182
void Wavefunction::init(void)
183
{
Francois Gygi committed
184
  // initialize all SlaterDets with lowest energy plane waves
185 186 187 188
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
189
      sd_[ispin][ikp]->init();
190 191
    }
  }
Francois Gygi committed
192
}
193

Francois Gygi committed
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::clear(void)
{
  // initialize all SlaterDets with zero
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
      sd_[ispin][ikp]->c().clear();
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::reset(void)
{
  // reset to single kpoint, ecut=0
  deallocate();
  nel_ = 0;
  nempty_ = 0;
  nspin_ = 1;
  deltaspin_ = 0;
  kpoint_.resize(1);
  kpoint_[0] = D3vector(0,0,0);
  weight_.resize(1);
  weight_[0] = 1.0;
  compute_nst();
  allocate();
}

Francois Gygi committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
////////////////////////////////////////////////////////////////////////////////
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)
{
245
  if ( nel == nel_ ) return;
Francois Gygi committed
246 247 248 249 250
  if ( nel < 0 )
  {
    cout << " Wavefunction::set_nel: nel < 0" << endl;
    return;
  }
251

Francois Gygi committed
252 253
  nel_ = nel;
  compute_nst();
254
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
255
}
256

Francois Gygi committed
257 258 259
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nempty(int nempty)
{
260
  if ( nempty == nempty_ ) return;
Francois Gygi committed
261 262 263 264 265 266 267
  if ( nempty < 0 )
  {
    cout << " Wavefunction::set_nempty: negative value" << endl;
    return;
  }
  nempty_ = nempty;
  compute_nst();
268
  update_occ(0.0);
269
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
270 271 272 273 274 275 276
}

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

Francois Gygi committed
278 279
  deallocate();
  nspin_ = nspin;
280

281
  compute_nst();
Francois Gygi committed
282
  allocate();
283
  resize(cell_,refcell_,ecut_);
284
  init();
285
  update_occ(0.0);
286 287 288 289 290 291 292
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_deltaspin(int deltaspin)
{
  if ( deltaspin == deltaspin_ ) return;

293 294 295 296 297 298 299 300 301 302 303
  // check if value of deltaspin would result in nst[1] < 0
  // nst_[1] = nel_ / 2 - deltaspin_ + nempty_;
  if ( nel_ / 2 - deltaspin + nempty_ < 0 )
  {
    if ( ctxt_.onpe0() )
    {
      cout << " Wavefunction::set_deltaspin: nel+nempty too small" << endl;
      cout << " Wavefunction::set_deltaspin: cannot set deltaspin" << endl;
      return;
    }
  }
304 305 306 307
  deallocate();
  deltaspin_ = deltaspin;
  compute_nst();
  allocate();
308
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
309
  init();
310
  update_occ(0.0);
Francois Gygi committed
311 312
}

313 314 315 316 317 318 319 320 321 322 323
////////////////////////////////////////////////////////////////////////////////
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;

324
  spincontext_ = new Context(ctxt_.comm(),npr,npc);
325 326 327 328
  kpcontext_ = new Context(*spincontext_);
  sdcontext_ = new Context(*kpcontext_);
}

Francois Gygi committed
329 330 331 332 333
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::set_nrowmax(int n)
{
  if ( n > ctxt_.size() )
  {
334 335
    if ( ctxt_.onpe0() )
      cout << " Wavefunction::set_nrowmax: nrowmax > ctxt_.size()" << endl;
Francois Gygi committed
336 337
    return;
  }
338

Francois Gygi committed
339
  deallocate();
340 341 342
  delete spincontext_;
  delete kpcontext_;
  delete sdcontext_;
Francois Gygi committed
343
  nrowmax_ = n;
344
  create_contexts();
345
  compute_nst();
Francois Gygi committed
346
  allocate();
347
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
348
  init();
Francois Gygi committed
349
}
350

Francois Gygi committed
351 352 353 354 355
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::add_kpoint(D3vector kpoint, double weight)
{
  for ( int i = 0; i < kpoint_.size(); i++ )
  {
Francois Gygi committed
356
    if ( length(kpoint - kpoint_[i]) < 1.e-6 )
Francois Gygi committed
357
    {
Francois Gygi committed
358 359
      if ( ctxt_.onpe0() )
        cout << " Wavefunction::add_kpoint: kpoint already defined"
Francois Gygi committed
360
           << endl;
Francois Gygi committed
361
      return;
Francois Gygi committed
362 363
    }
  }
364

Francois Gygi committed
365 366
  kpoint_.push_back(kpoint);
  weight_.push_back(weight);
367

368 369
  const int nkp = kpoint_.size();
  sd_.resize(nspin_);
370 371 372 373 374
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    sd_[ispin].push_back(new SlaterDet(*sdcontext_,kpoint_[nkp-1]));
    sd_[ispin][nkp-1]->resize(cell_,refcell_,ecut_,nst_[ispin]);
  }
375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391

  if ( nspin_ == 1 )
  {
    sd_[0][nkp-1]->update_occ(nel_,nspin_);
  }
  else if ( nspin_ == 2 )
  {
    const int nocc_up = (nel_+1)/2+deltaspin_;
    const int nocc_dn = nel_/2 - deltaspin_;
    sd_[0][nkp-1]->update_occ(nocc_up,nspin_);
    sd_[1][nkp-1]->update_occ(nocc_dn,nspin_);
  }
  else
  {
    // incorrect value of nspin_
    assert(false);
  }
Francois Gygi committed
392 393 394 395 396
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::del_kpoint(D3vector kpoint)
{
Francois Gygi committed
397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
  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_);
Francois Gygi committed
424
  init();
Francois Gygi committed
425
  update_occ(0.0);
Francois Gygi committed
426 427 428 429 430 431 432 433 434
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::randomize(double amplitude)
{
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
435
      sd_[ispin][ikp]->randomize(amplitude);
Francois Gygi committed
436 437 438 439 440 441 442
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::update_occ(double temp)
{
Francois Gygi committed
443 444
  // update occupation numbers using eigenvalues in SlaterDet
  if ( temp == 0.0 )
Francois Gygi committed
445
  {
Francois Gygi committed
446 447
    // zero temperature
    if ( nspin_ == 1 )
Francois Gygi committed
448
    {
Francois Gygi committed
449
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
450
      {
Francois Gygi committed
451
        sd_[0][ikp]->update_occ(nel_,nspin_);
Francois Gygi committed
452 453
      }
    }
Francois Gygi committed
454
    else if ( nspin_ == 2 )
Francois Gygi committed
455
    {
Francois Gygi committed
456 457 458 459
      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
460 461
        sd_[0][ikp]->update_occ(nocc_up,nspin_);
        sd_[1][ikp]->update_occ(nocc_dn,nspin_);
Francois Gygi committed
462
      }
Francois Gygi committed
463 464 465
    }
    else
    {
Francois Gygi committed
466 467
      // incorrect value of nspin_
      assert(false);
Francois Gygi committed
468
    }
Francois Gygi committed
469 470 471 472 473
  }
  else
  {
    // finite temperature
    const double eVolt = 0.036749023; // 1 eV in Hartree
474
    const int maxiter = 500;
475

Francois Gygi committed
476 477
    // loop to find value of mu
    double mu = 0.0;
478
    double dmu = 2.0 * eVolt;
Francois Gygi committed
479 480 481 482 483
    const double totalcharge = (double) nel_;
    enum direction { up, down };
    direction dir = up;

    double rhosum = 0.0;
Francois Gygi committed
484 485 486 487
    for ( int ispin = 0; ispin < nspin_; ispin++ )
    {
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
      {
Francois Gygi committed
488
        sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
489
        rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
490 491
      }
    }
492

Francois Gygi committed
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
    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
514
          sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
515
          rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
516 517 518
        }
      }
    }
519

Francois Gygi committed
520 521 522 523 524 525
    if ( niter == maxiter )
    {
      cout << "Wavefunction::update_occ: mu did not converge in "
           << maxiter << " iterations" << endl;
      ctxt_.abort(1);
    }
Francois Gygi committed
526

Francois Gygi committed
527
    if ( ctxt_.onpe0() )
Francois Gygi committed
528
    {
529 530 531 532
      cout << " Wavefunction::update_occ: sum = "
           << rhosum << endl;
      cout << " Wavefunction::update_occ: mu = "
           << setprecision(4) << mu / eVolt << " eV" << endl;
Francois Gygi committed
533 534 535

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

537
      cout << " Wavefunction::update_occ: occupation numbers" << endl;
Francois Gygi committed
538
      for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
539
      {
Francois Gygi committed
540
        for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
541
        {
542
          cout << " k = " << kpoint_[ikp] << endl;
Francois Gygi committed
543
          for ( int n = 0; n < sd_[ispin][ikp]->nst(); n++ )
Francois Gygi committed
544
          {
Francois Gygi committed
545 546
            cout << setw(7) << setprecision(4) << sd_[ispin][ikp]->occ(n);
            if ( ( n%10 ) == 9 ) cout << endl;
Francois Gygi committed
547
          }
548 549
          if ( sd_[ispin][ikp]->nst() % 10 != 0 )
            cout << endl;
Francois Gygi committed
550 551 552 553 554 555 556
        }
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
557
void Wavefunction::gram(void)
Francois Gygi committed
558
{
Francois Gygi committed
559
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
560 561 562
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
563
      sd_[ispin][ikp]->gram();
Francois Gygi committed
564 565
    }
  }
Francois Gygi committed
566 567 568 569 570 571 572
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::riccati(Wavefunction& wf)
{
  assert(wf.context() == ctxt_);
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
573 574 575
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
576
      sd_[ispin][ikp]->riccati(*wf.sd_[ispin][ikp]);
Francois Gygi committed
577 578 579 580 581
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
582
void Wavefunction::align(Wavefunction& wf)
Francois Gygi committed
583
{
Francois Gygi committed
584
  assert(wf.context() == ctxt_);
Francois Gygi committed
585 586 587 588
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
589
      sd_[ispin][ikp]->align(*wf.sd_[ispin][ikp]);
Francois Gygi committed
590 591 592 593 594
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
595
complex<double> Wavefunction::dot(const Wavefunction& wf) const
Francois Gygi committed
596 597
{
  assert(wf.context() == ctxt_);
Francois Gygi committed
598
  complex<double> sum = 0.0;
Francois Gygi committed
599 600 601 602
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
603
      sum += sd_[ispin][ikp]->dot(*wf.sd_[ispin][ikp]);
Francois Gygi committed
604 605 606 607 608 609 610 611 612
    }
  }
  return sum;
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::diag(Wavefunction& dwf, bool eigvec)
{
  // subspace diagonalization of <*this | dwf>
613
  // if eigvec==true, eigenvectors are computed and stored in *this, dwf is
Francois Gygi committed
614 615 616
  // overwritten
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
617
    if ( nst_[ispin] > 0 )
Francois Gygi committed
618
    {
619
      for ( int ikp = 0; ikp < nkp(); ikp++ )
Francois Gygi committed
620
      {
621 622 623 624 625 626
        // compute eigenvalues
        if ( sd(ispin,ikp)->basis().real() )
        {
          // proxy real matrices c, cp
          DoubleMatrix c(sd(ispin,ikp)->c());
          DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
627

628
          DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
629

630 631 632 633
          // 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);
634

635 636 637
          // cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
          //      << endl;
          // cout << h;
638

639
#if 1
640 641 642 643 644 645 646 647 648 649 650 651 652 653 654
          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);
          }
#else
          vector<double> w(h.m());
Francois Gygi committed
655
          DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
656 657 658 659 660 661 662 663 664 665
          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);
          }
#endif
          // set eigenvalues in SlaterDet
          sd(ispin,ikp)->set_eig(w);
Francois Gygi committed
666 667 668
        }
        else
        {
669 670 671 672 673 674 675 676 677 678
          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);
          // cout << " Hamiltonian at k = "
          //      << sd(ispin,ikp)->kpoint() << endl;
          // cout << h;
          valarray<double> w(h.m());
          if ( eigvec )
          {
679
#if DEBUG
680
            ComplexMatrix hcopy(h);
681
#endif
682 683 684 685
            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);
686
#if DEBUG
687 688 689 690 691 692 693 694 695 696 697 698 699 700
            // 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
701
#endif
702 703 704 705 706 707 708
          }
          else
          {
            h.heev('l',w);
          }
          // set eigenvalues in SlaterDet
          sd(ispin,ikp)->set_eig(w);
709
        }
Francois Gygi committed
710 711 712 713 714 715
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
716
void Wavefunction::print(ostream& os, string encoding, string tag) const
Francois Gygi committed
717
{
718 719 720 721 722 723
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
724
    os << setprecision(10);
725 726 727 728 729 730 731 732
    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;
  }
733

Francois Gygi committed
734 735 736
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
737
      sd_[ispin][ikp]->print(os,encoding,weight_[ikp],ispin,nspin_);
Francois Gygi committed
738
  }
739

740 741
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
742 743 744
}

////////////////////////////////////////////////////////////////////////////////
745
void Wavefunction::write(SharedFilePtr& sfp, string encoding, string tag) const
Francois Gygi committed
746
{
747
  sfp.sync();
748

Francois Gygi committed
749 750 751 752 753 754 755
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
756
    os << setprecision(10);
Francois Gygi committed
757 758 759 760
    os << "<domain a=\""
       << cell_.a(0) << "\"\n        b=\""
       << cell_.a(1) << "\"\n        c=\""
       << cell_.a(2) << "\"/>" << endl;
761 762 763 764 765 766 767
    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
768 769 770 771
    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());
772
    int len = str.size();
773
#if USE_MPI
774
    MPI_Status status;
775 776
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
777
    if ( err != 0 )
778
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
779
    sfp.advance(len);
780 781 782
#else
    sfp.file().write(str.c_str(),len);
#endif
Francois Gygi committed
783 784 785 786 787
  }

  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
788
    {
789
      sd_[ispin][ikp]->write(sfp,encoding,weight_[ikp],ispin,nspin_);
790
    }
Francois Gygi committed
791
  }
792

793
  sfp.sync();
794

Francois Gygi committed
795 796 797 798 799
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "</" << tag << ">" << endl;
    string str(os.str());
800
    int len = str.size();
801
#if USE_MPI
802
    MPI_Status status;
803 804
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
805
    if ( err != 0 )
806
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
807
    sfp.advance(len);
808 809 810
#else
    sfp.file().write(str.c_str(),len);
#endif
Francois Gygi committed
811 812 813 814
  }
}

////////////////////////////////////////////////////////////////////////////////
815
void Wavefunction::info(ostream& os, string tag) const
Francois Gygi committed
816
{
817 818 819 820 821 822
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
823
    os.setf(ios::fixed,ios::floatfield);
824
    os << "<cell a=\""
825 826
       << setprecision(6) << cell_.a(0) << "\"\n      b=\""
       << cell_.a(1) << "\"\n      c=\""
827
       << cell_.a(2) << "\"/>" << endl;
828
    os << " reciprocal lattice vectors" << endl
Francois Gygi committed
829 830 831
       << setprecision(6)
       << " " << cell_.b(0) << endl
       << " " << cell_.b(1) << endl
832
       << " " << cell_.b(2) << endl;
833
    os << "<refcell a=\""
834 835
       << refcell_.a(0) << "\"\n         b=\""
       << refcell_.a(1) << "\"\n         c=\""
836
       << refcell_.a(2) << "\"/>" << endl;
837 838 839 840
    os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
       <<      " ny=\"" << sd_[0][0]->basis().np(1) << "\""
       <<      " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
  }
841

Francois Gygi committed
842 843 844
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
845
    {
846
      if ( ctxt_.onpe0() )
847 848
        cout << " kpoint: " << kpoint_[ikp] << " weight: " << weight_[ikp]
             << endl;
Francois Gygi committed
849
      sd_[ispin][ikp]->info(os);
Francois Gygi committed
850
    }
Francois Gygi committed
851
  }
852

853 854
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
855 856 857
}

////////////////////////////////////////////////////////////////////////////////
858
ostream& operator<<(ostream& os, const Wavefunction& wf)
Francois Gygi committed
859 860 861 862 863 864 865 866 867 868 869 870 871
{
  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_);
872
  assert(nrowmax_ == wf.nrowmax_);
Francois Gygi committed
873 874 875
  assert(deltaspin_ == wf.deltaspin_);
  assert(refcell_ == wf.refcell_);
  assert(ecut_ == wf.ecut_);
876

877
  cell_ = wf.cell_;
Francois Gygi committed
878 879
  weight_ = wf.weight_;
  kpoint_ = wf.kpoint_;
880

881 882 883 884
  assert(*spincontext_ == *wf.spincontext_);
  assert(*kpcontext_   == *wf.kpcontext_);
  assert(*sdcontext_   == *wf.sdcontext_);

Francois Gygi committed
885 886 887 888
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
889
      *sd_[ispin][ikp] = *wf.sd_[ispin][ikp];
Francois Gygi committed
890 891 892 893
    }
  }
  return *this;
}