Wavefunction.C 26.7 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 "FourierTransform.h"
22
#include "jacobi.h"
23
#include "SharedFilePtr.h"
Francois Gygi committed
24 25 26 27 28 29
#include <vector>
#include <iomanip>
#include <sstream>
using namespace std;

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

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

Francois Gygi committed
52
  compute_nst();
53 54

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

60
  allocate();
61

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

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

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

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

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

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

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

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

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

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

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

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

180
}
181

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

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

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

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

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

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

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

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

294 295 296 297 298 299 300 301 302 303 304
  // 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;
    }
  }
305 306 307 308
  deallocate();
  deltaspin_ = deltaspin;
  compute_nst();
  allocate();
309
  resize(cell_,refcell_,ecut_);
Francois Gygi committed
310
  init();
311
  update_occ(0.0);
Francois Gygi committed
312 313
}

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

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

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

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

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

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

369 370
  const int nkp = kpoint_.size();
  sd_.resize(nspin_);
371 372 373 374 375
  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]);
  }
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392

  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
393 394 395 396 397
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::del_kpoint(D3vector kpoint)
{
Francois Gygi committed
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 424
  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
425
  init();
Francois Gygi committed
426
  update_occ(0.0);
Francois Gygi committed
427 428
}

429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::move_kpoint(D3vector kpoint, D3vector newkpoint)
{
  bool found = false;
  int ikp = 0;
  while ( !found && ikp < kpoint_.size() )
  {
    if ( length(kpoint_[ikp] - kpoint) < 1.e-6 )
    {
      found = true;
    }
    else
    {
      ikp++;
    }
  }
  if ( !found )
  {
    if ( ctxt_.onpe0() )
      cout << " Wavefunction::move_kpoint: no such kpoint"
         << endl;
    return;
  }
  // check if newkpoint coincides with existing kpoint
  for ( int i = 0; i < kpoint_.size(); i++ )
  {
    if ( length(newkpoint - kpoint_[i]) < 1.e-6 )
    {
      if ( ctxt_.onpe0() )
        cout << " Wavefunction::move_kpoint: kpoint already defined "
             << "at newkpoint position"
             << endl;
      return;
    }
  }

  // copy wavefunctions from old SlaterDet at kpoint to new SlaterDet
  // at newkpoint
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    // create new SlaterDet at newkpoint
    SlaterDet *sd = sd_[ispin][ikp];
    SlaterDet *sdn = new SlaterDet(*sdcontext_,newkpoint);
    sdn->resize(cell_,refcell_,ecut_,nst_[ispin]);
    sdn->init();
    // copy wave functions from old to new SlaterDet
    const Basis& basis = sd_[ispin][ikp]->basis();
    const Basis& newbasis = sdn->basis();
    // transform all states to real space and interpolate
    int np0 = max(basis.np(0),newbasis.np(0));
    int np1 = max(basis.np(1),newbasis.np(1));
    int np2 = max(basis.np(2),newbasis.np(2));
    FourierTransform ft1(basis,np0,np1,np2);
    FourierTransform ft2(newbasis,np0,np1,np2);
    // allocate real-space grid
    valarray<complex<double> > tmpr(ft1.np012loc());
    for ( int n = 0; n < sd->nstloc(); n++ )
    {
      for ( int i = 0; i < tmpr.size(); i++ )
        tmpr[i] = 0.0;
      ComplexMatrix& c = sd->c();
      ComplexMatrix& cn = sdn->c();
      ft1.backward(c.valptr(n*c.mloc()),&tmpr[0]);
      // if the new kpoint is Gamma, remove the phase of the wavefunction
      // using the phase of the first element of the array
      if ( newbasis.real() )
      {
        double arg0 = arg(tmpr[0]);
        // broadcast argument found on task 0
        MPI_Bcast(&arg0,1,MPI_DOUBLE,0,basis.comm());
        complex<double> z = polar(1.0,-arg0);
        for ( int i = 0; i < tmpr.size(); i++ )
          tmpr[i] *= z;
      }
      ft2.forward(&tmpr[0], cn.valptr(n*cn.mloc()));
    }

    // delete old SlaterDet
    delete sd_[ispin][ikp];
    // reassign pointer
    sd_[ispin][ikp] = sdn;
  }
  kpoint_[ikp] = newkpoint;

  if ( nspin_ == 1 )
  {
    sd_[0][ikp]->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][ikp]->update_occ(nocc_up,nspin_);
    sd_[1][ikp]->update_occ(nocc_dn,nspin_);
  }
  else
  {
    // incorrect value of nspin_
    assert(false);
  }
}

Francois Gygi committed
531 532 533 534 535 536 537
////////////////////////////////////////////////////////////////////////////////
void Wavefunction::randomize(double amplitude)
{
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
538
      sd_[ispin][ikp]->randomize(amplitude);
Francois Gygi committed
539 540 541 542 543 544 545
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::update_occ(double temp)
{
Francois Gygi committed
546 547
  // update occupation numbers using eigenvalues in SlaterDet
  if ( temp == 0.0 )
Francois Gygi committed
548
  {
Francois Gygi committed
549 550
    // zero temperature
    if ( nspin_ == 1 )
Francois Gygi committed
551
    {
Francois Gygi committed
552
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
553
      {
Francois Gygi committed
554
        sd_[0][ikp]->update_occ(nel_,nspin_);
Francois Gygi committed
555 556
      }
    }
Francois Gygi committed
557
    else if ( nspin_ == 2 )
Francois Gygi committed
558
    {
Francois Gygi committed
559 560 561 562
      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
563 564
        sd_[0][ikp]->update_occ(nocc_up,nspin_);
        sd_[1][ikp]->update_occ(nocc_dn,nspin_);
Francois Gygi committed
565
      }
Francois Gygi committed
566 567 568
    }
    else
    {
Francois Gygi committed
569 570
      // incorrect value of nspin_
      assert(false);
Francois Gygi committed
571
    }
Francois Gygi committed
572 573 574 575 576
  }
  else
  {
    // finite temperature
    const double eVolt = 0.036749023; // 1 eV in Hartree
577
    const int maxiter = 500;
578

Francois Gygi committed
579 580
    // loop to find value of mu
    double mu = 0.0;
581
    double dmu = 2.0 * eVolt;
Francois Gygi committed
582 583 584 585 586
    const double totalcharge = (double) nel_;
    enum direction { up, down };
    direction dir = up;

    double rhosum = 0.0;
Francois Gygi committed
587 588 589 590
    for ( int ispin = 0; ispin < nspin_; ispin++ )
    {
      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
      {
Francois Gygi committed
591
        sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
592
        rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
593 594
      }
    }
595

Francois Gygi committed
596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616
    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
617
          sd_[ispin][ikp]->update_occ(nspin_,mu,temp);
618
          rhosum += weight_[ikp] * sd_[ispin][ikp]->total_charge();
Francois Gygi committed
619 620 621
        }
      }
    }
622

Francois Gygi committed
623 624 625 626 627 628
    if ( niter == maxiter )
    {
      cout << "Wavefunction::update_occ: mu did not converge in "
           << maxiter << " iterations" << endl;
      ctxt_.abort(1);
    }
Francois Gygi committed
629

Francois Gygi committed
630
    if ( ctxt_.onpe0() )
Francois Gygi committed
631
    {
632 633 634 635
      cout << " Wavefunction::update_occ: sum = "
           << rhosum << endl;
      cout << " Wavefunction::update_occ: mu = "
           << setprecision(4) << mu / eVolt << " eV" << endl;
Francois Gygi committed
636 637 638

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

640
      cout << " Wavefunction::update_occ: occupation numbers" << endl;
Francois Gygi committed
641
      for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
642
      {
Francois Gygi committed
643
        for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
644
        {
645
          cout << " k = " << kpoint_[ikp] << endl;
Francois Gygi committed
646
          for ( int n = 0; n < sd_[ispin][ikp]->nst(); n++ )
Francois Gygi committed
647
          {
Francois Gygi committed
648 649
            cout << setw(7) << setprecision(4) << sd_[ispin][ikp]->occ(n);
            if ( ( n%10 ) == 9 ) cout << endl;
Francois Gygi committed
650
          }
651 652
          if ( sd_[ispin][ikp]->nst() % 10 != 0 )
            cout << endl;
Francois Gygi committed
653 654 655 656 657 658 659
        }
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
660
void Wavefunction::gram(void)
Francois Gygi committed
661
{
Francois Gygi committed
662
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
663 664 665
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
666
      sd_[ispin][ikp]->gram();
Francois Gygi committed
667 668
    }
  }
Francois Gygi committed
669 670 671 672 673 674 675
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::riccati(Wavefunction& wf)
{
  assert(wf.context() == ctxt_);
  for ( int ispin = 0; ispin < nspin_; ispin++ )
Francois Gygi committed
676 677 678
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
679
      sd_[ispin][ikp]->riccati(*wf.sd_[ispin][ikp]);
Francois Gygi committed
680 681 682 683 684
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
685
void Wavefunction::align(Wavefunction& wf)
Francois Gygi committed
686
{
Francois Gygi committed
687
  assert(wf.context() == ctxt_);
Francois Gygi committed
688 689 690 691
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
692
      sd_[ispin][ikp]->align(*wf.sd_[ispin][ikp]);
Francois Gygi committed
693 694 695 696 697
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
698
complex<double> Wavefunction::dot(const Wavefunction& wf) const
Francois Gygi committed
699 700
{
  assert(wf.context() == ctxt_);
Francois Gygi committed
701
  complex<double> sum = 0.0;
Francois Gygi committed
702 703 704 705
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
706
      sum += weight_[ikp] * sd_[ispin][ikp]->dot(*wf.sd_[ispin][ikp]);
Francois Gygi committed
707 708 709 710 711 712 713 714 715
    }
  }
  return sum;
}

////////////////////////////////////////////////////////////////////////////////
void Wavefunction::diag(Wavefunction& dwf, bool eigvec)
{
  // subspace diagonalization of <*this | dwf>
716
  // if eigvec==true, eigenvectors are computed and stored in *this, dwf is
Francois Gygi committed
717 718 719
  // overwritten
  for ( int ispin = 0; ispin < nspin(); ispin++ )
  {
720
    if ( nst_[ispin] > 0 )
Francois Gygi committed
721
    {
722
      for ( int ikp = 0; ikp < nkp(); ikp++ )
Francois Gygi committed
723
      {
724 725 726 727 728 729
        // 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());
730

731
          DoubleMatrix h(c.context(),c.n(),c.n(),c.nb(),c.nb());
732

733 734 735 736
          // 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);
737

738 739 740
          // cout << " Hamiltonian at k = " << sd(ispin,ikp)->kpoint()
          //      << endl;
          // cout << h;
741

742
#if 1
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757
          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
758
          DoubleMatrix z(c.context(),c.n(),c.n(),c.nb(),c.nb());
759 760 761 762 763 764 765 766 767 768
          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
769 770 771
        }
        else
        {
772 773 774 775 776 777 778 779 780 781
          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 )
          {
782
#if DEBUG
783
            ComplexMatrix hcopy(h);
784
#endif
785 786 787 788
            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);
789
#if DEBUG
790 791 792 793 794 795 796 797 798 799 800 801 802 803
            // 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
804
#endif
805 806 807 808 809 810 811
          }
          else
          {
            h.heev('l',w);
          }
          // set eigenvalues in SlaterDet
          sd(ispin,ikp)->set_eig(w);
812
        }
Francois Gygi committed
813 814 815 816 817 818
      }
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
819
void Wavefunction::print(ostream& os, string encoding, string tag) const
Francois Gygi committed
820
{
821 822 823 824 825 826
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
827
    os << setprecision(10);
828 829 830 831 832 833 834 835
    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;
  }
836

Francois Gygi committed
837 838 839
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
840
      sd_[ispin][ikp]->print(os,encoding,weight_[ikp],ispin,nspin_);
Francois Gygi committed
841
  }
842

843 844
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
845 846 847
}

////////////////////////////////////////////////////////////////////////////////
848
void Wavefunction::write(SharedFilePtr& sfp, string encoding, string tag) const
Francois Gygi committed
849
{
850
  sfp.sync();
851

Francois Gygi committed
852 853 854 855 856 857 858
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
859
    os << setprecision(10);
Francois Gygi committed
860 861 862 863
    os << "<domain a=\""
       << cell_.a(0) << "\"\n        b=\""
       << cell_.a(1) << "\"\n        c=\""
       << cell_.a(2) << "\"/>" << endl;
864 865 866 867 868 869 870
    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
871 872 873 874
    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());
875
    int len = str.size();
876
#if USE_MPI
877
    MPI_Status status;
878 879
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
880
    if ( err != 0 )
881
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
882
    sfp.advance(len);
883 884 885
#else
    sfp.file().write(str.c_str(),len);
#endif
Francois Gygi committed
886 887 888 889 890
  }

  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
891
    {
892
      sd_[ispin][ikp]->write(sfp,encoding,weight_[ikp],ispin,nspin_);
893
    }
Francois Gygi committed
894
  }
895

896
  sfp.sync();
897

Francois Gygi committed
898 899 900 901 902
  if ( ctxt_.onpe0() )
  {
    ostringstream os;
    os << "</" << tag << ">" << endl;
    string str(os.str());
903
    int len = str.size();
904
#if USE_MPI
905
    MPI_Status status;
906 907
    int err = MPI_File_write_at(sfp.file(),sfp.mpi_offset(),(void*)str.c_str(),
              len,MPI_CHAR,&status);
908
    if ( err != 0 )
909
      cout << " Wavefunction::write: error in MPI_File_write" << endl;
910
    sfp.advance(len);
911 912 913
#else
    sfp.file().write(str.c_str(),len);
#endif
Francois Gygi committed
914 915 916 917
  }
}

////////////////////////////////////////////////////////////////////////////////
918
void Wavefunction::info(ostream& os, string tag) const
Francois Gygi committed
919
{
920 921 922 923 924 925
  if ( ctxt_.onpe0() )
  {
    os << "<" << tag << " ecut=\"" << ecut_ << "\""
       << " nspin=\"" << nspin_ << "\""
       << " nel=\"" << nel_ << "\""
       << " nempty=\"" << nempty_ << "\">" << endl;
926
    os.setf(ios::fixed,ios::floatfield);
927
    os << "<cell a=\""
928 929
       << setprecision(6) << cell_.a(0) << "\"\n      b=\""
       << cell_.a(1) << "\"\n      c=\""
930
       << cell_.a(2) << "\"/>" << endl;
931
    os << " reciprocal lattice vectors" << endl
Francois Gygi committed
932 933 934
       << setprecision(6)
       << " " << cell_.b(0) << endl
       << " " << cell_.b(1) << endl
935
       << " " << cell_.b(2) << endl;
936
    os << "<refcell a=\""
937 938
       << refcell_.a(0) << "\"\n         b=\""
       << refcell_.a(1) << "\"\n         c=\""
939
       << refcell_.a(2) << "\"/>" << endl;
940 941 942 943
    os << "<grid nx=\"" << sd_[0][0]->basis().np(0) << "\""
       <<      " ny=\"" << sd_[0][0]->basis().np(1) << "\""
       <<      " nz=\"" << sd_[0][0]->basis().np(2) << "\"/>" << endl;
  }
944

Francois Gygi committed
945 946 947
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
Francois Gygi committed
948
    {
949
      if ( ctxt_.onpe0() )
950 951
        cout << " kpoint: " << kpoint_[ikp] << " weight: " << weight_[ikp]
             << endl;
Francois Gygi committed
952
      sd_[ispin][ikp]->info(os);
Francois Gygi committed
953
    }
Francois Gygi committed
954
  }
955

956 957
  if ( ctxt_.onpe0() )
    os << "</" << tag << ">" << endl;
Francois Gygi committed
958 959 960
}

////////////////////////////////////////////////////////////////////////////////
961
ostream& operator<<(ostream& os, const Wavefunction& wf)
Francois Gygi committed
962 963 964 965 966 967 968 969 970 971 972 973 974
{
  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_);
975
  assert(nrowmax_ == wf.nrowmax_);
Francois Gygi committed
976 977 978
  assert(deltaspin_ == wf.deltaspin_);
  assert(refcell_ == wf.refcell_);
  assert(ecut_ == wf.ecut_);
979

980
  cell_ = wf.cell_;
Francois Gygi committed
981 982
  weight_ = wf.weight_;
  kpoint_ = wf.kpoint_;
983

984 985 986 987
  assert(*spincontext_ == *wf.spincontext_);
  assert(*kpcontext_   == *wf.kpcontext_);
  assert(*sdcontext_   == *wf.sdcontext_);

Francois Gygi committed
988 989 990 991
  for ( int ispin = 0; ispin < nspin_; ispin++ )
  {
    for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
    {
Francois Gygi committed
992
      *sd_[ispin][ikp] = *wf.sd_[ispin][ikp];
Francois Gygi committed
993 994 995 996
    }
  }
  return *this;
}