Basis.C 21.9 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
// Basis.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18
// $Id: Basis.C,v 1.22 2008-09-08 15:56:18 fgygi Exp $
Francois Gygi committed
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

#include "Basis.h"
#include "Context.h"

#include <cmath>
#include <cassert>
#include <cstdlib>
#include <vector>
#include <set>
#include <algorithm>
#include <functional>
#include <iostream>
#include <iomanip>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
double Basis::localmemsize(void) const
36 37
{
  return
Francois Gygi committed
38 39
  5.0 * (nprow_*nrods_*sizeof(int)) // x[ipe][irod]
  + localsize_[myrow_] * (3.0*sizeof(int) + 10 * sizeof(double));
Francois Gygi committed
40
}
Francois Gygi committed
41
double Basis::memsize(void) const { return nprow_*localmemsize(); }
Francois Gygi committed
42

Francois Gygi committed
43
const Context& Basis::context(void) const { return ctxt_; }
Francois Gygi committed
44

Francois Gygi committed
45 46 47 48 49
const UnitCell& Basis::cell() const { return cell_; }
const UnitCell& Basis::refcell() const { return refcell_; }
int Basis::idxmin(int i) const { return idxmin_[i]; }
int Basis::idxmax(int i) const { return idxmax_[i]; }
double Basis::ecut() const { return ecut_; }
Francois Gygi committed
50

Francois Gygi committed
51 52 53 54 55
int Basis::size() const { return size_; }
int Basis::localsize() const { return localsize_[myrow_]; }
int Basis::localsize(int ipe) const { return localsize_[ipe]; }
int Basis::maxlocalsize() const { return maxlocalsize_; }
int Basis::minlocalsize() const { return minlocalsize_; }
Francois Gygi committed
56

Francois Gygi committed
57 58 59
int Basis::nrods() const { return nrods_; }
int Basis::nrod_loc() const { return nrod_loc_[myrow_]; }
int Basis::nrod_loc(int ipe) const { return nrod_loc_[ipe]; }
Francois Gygi committed
60 61

int Basis::rod_h(int irod) const
Francois Gygi committed
62
{ return rod_h_[myrow_][irod]; }
Francois Gygi committed
63
int Basis::rod_h(int ipe, int irod) const
Francois Gygi committed
64
{ return rod_h_[ipe][irod]; }
Francois Gygi committed
65

Francois Gygi committed
66 67
int Basis::rod_k(int irod) const { return rod_k_[myrow_][irod]; }
int Basis::rod_k(int ipe, int irod) const { return rod_k_[ipe][irod]; }
Francois Gygi committed
68

Francois Gygi committed
69 70
int Basis::rod_lmin(int irod) const { return rod_lmin_[myrow_][irod]; }
int Basis::rod_lmin(int ipe, int irod) const { return rod_lmin_[ipe][irod]; }
Francois Gygi committed
71 72

// size of rod irod on current process
Francois Gygi committed
73 74
int Basis::rod_size(int irod) const { return rod_size_[myrow_][irod]; }
int Basis::rod_size(int ipe, int irod) const { return rod_size_[ipe][irod]; }
Francois Gygi committed
75 76

// position of first elem. of rod irod in the local list of g vectors
Francois Gygi committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
int Basis::rod_first(int irod) const { return rod_first_[myrow_][irod]; }
int Basis::rod_first(int ipe, int irod) const { return rod_first_[ipe][irod]; }

int    Basis::idx(int i) const   { return idx_[i]; }
double Basis::g(int i) const     { return g_[i]; }
double Basis::kpg(int i) const   { return kpg_[i]; }
double Basis::gi(int i) const    { return gi_[i]; }
double Basis::kpgi(int i) const  { return kpgi_[i]; }
double Basis::g2(int i) const    { return g2_[i]; }
double Basis::kpg2(int i) const  { return kpg2_[i]; }
double Basis::g2i(int i) const   { return g2i_[i]; }
double Basis::kpg2i(int i) const { return kpg2i_[i]; }
double Basis::gx(int i) const    { return gx_[i]; }
double Basis::kpgx(int i) const  { return kpgx_[i]; }
int    Basis::isort(int i) const { return isort_loc[i]; }

const int*    Basis::idx_ptr(void) const   { return &(idx_[0]); }
const double* Basis::g_ptr(void)  const    { return &(g_[0]); }
const double* Basis::kpg_ptr(void)  const  { return &(kpg_[0]); }
const double* Basis::gi_ptr(void) const    { return &(gi_[0]); }
const double* Basis::kpgi_ptr(void) const  { return &(kpgi_[0]); }
const double* Basis::g2_ptr(void) const    { return &(g2_[0]); }
const double* Basis::kpg2_ptr(void) const  { return &(kpg2_[0]); }
const double* Basis::g2i_ptr(void) const   { return &(g2i_[0]); }
const double* Basis::kpg2i_ptr(void) const { return &(kpg2i_[0]); }
102
const double* Basis::gx_ptr(int j) const
Francois Gygi committed
103
{ return &(gx_[j*localsize_[myrow_]]); }
Francois Gygi committed
104
const double* Basis::kpgx_ptr(int j) const
Francois Gygi committed
105
{ return &(kpgx_[j*localsize_[myrow_]]); }
Francois Gygi committed
106 107

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
108
bool Basis::factorizable(int n) const
Francois Gygi committed
109 110
{
  // next lines: use AIX criterion for all platforms (AIX and fftw)
111

Francois Gygi committed
112 113 114
//#if AIX

  // Acceptable lengths for FFTs in the ESSL library:
115 116
  // n = (2^h) (3^i) (5^j) (7^k) (11^m) for n <= 37748736
  // where:
Francois Gygi committed
117 118 119 120 121 122 123 124 125 126 127 128
  //  h = 1, 2, ..., 25
  //  i = 0, 1, 2
  //  j, k, m = 0, 1
  if ( n % 11 == 0 ) n /= 11;
  if ( n % 7 == 0 ) n /= 7;
  if ( n % 5 == 0 ) n /= 5;
  if ( n % 3 == 0 ) n /= 3;
  if ( n % 3 == 0 ) n /= 3;
  // the bound h <= 25 is not tested since 2^25 would cause
  // memory allocation problems
  while ( ( n % 2 == 0 ) ) n /= 2;
  return ( n == 1 );
129

Francois Gygi committed
130 131 132 133 134 135 136 137 138
// #else
//   while ( n % 5 == 0 ) n /= 5;
//   while ( n % 3 == 0 ) n /= 3;
//   while ( n % 2 == 0 ) n /= 2;
//   return ( n == 1 );
// #endif
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
139
int Basis::np(int i) const { return np_[i]; }
Francois Gygi committed
140 141

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
142
const D3vector Basis::kpoint(void) const { return kpoint_; }
Francois Gygi committed
143 144

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
145
bool Basis::real(void) const { return real_; }
Francois Gygi committed
146 147 148 149 150 151 152 153 154

inline double sqr( double x ) { return x*x; }
inline void swap(int &x, int &y) { int tmp = x; x = y; y = tmp; }

////////////////////////////////////////////////////////////////////////////////
class Rod
{
  // z-column of non-zero reciprocal lattice vectors
  int h_, k_, lmin_, size_;
155

Francois Gygi committed
156
  public:
157
  Rod(int h, int k, int lmin, int size) : h_(h), k_(k),
Francois Gygi committed
158
  lmin_(lmin), size_(size) {}
159

Francois Gygi committed
160 161 162 163 164 165 166 167 168 169 170 171 172 173
  int h(void) const { return h_; }
  int k(void) const { return k_; }
  int lmin(void) const { return lmin_; }
  int size(void) const { return size_; }
  bool operator< (const Rod& x) const
  {
    return size_ < x.size();
  }
};

////////////////////////////////////////////////////////////////////////////////
class Node
{
  int id_, nrods_, size_;
174

Francois Gygi committed
175 176 177
  public:
  Node() : id_(0), nrods_(0), size_(0) {}
  Node(int id) : id_(id), nrods_(0), size_(0) {}
178

Francois Gygi committed
179 180 181
  int id(void) const { return id_; }
  int nrods(void) const { return nrods_; }
  int size(void) const { return size_; }
182

Francois Gygi committed
183 184 185 186 187 188 189 190 191 192 193 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
  void addrod(const Rod& r)
  {
    nrods_++;
    size_ += r.size();
  }
  bool operator< (const Node& x) const
  {
    return size_ < x.size();
  }
};

////////////////////////////////////////////////////////////////////////////////
template <class T>
struct ptr_less
{
  public:
  bool operator() ( T* x, T* y ) { return *x < *y; }
};
template <class T>
struct ptr_greater
{
  public:
  bool operator() ( T* x, T* y ) { return *y < *x; }
};

////////////////////////////////////////////////////////////////////////////////
template <class T>
struct VectorLess
{
  // function object for indirect comparison of vector elements
  public:
  vector<T>& a_;
  VectorLess<T>(vector<T>& a) : a_(a) {};
  bool operator() (int i, int j) const
  {
    return a_[i] < a_[j];
  }
};
221

Francois Gygi committed
222
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
223
Basis::Basis(const Context& ctxt, D3vector kpoint) : ctxt_(ctxt)
Francois Gygi committed
224 225 226 227 228 229 230 231 232
{
  // Construct the default empty basis
  // cell and refcell are (0,0,0)
  myrow_ = ctxt.myrow();
  nprow_ = ctxt.nprow();

  ecut_ = 0.0;
  kpoint_ = kpoint;
  real_ = ( kpoint == D3vector(0.0,0.0,0.0) );
233

Francois Gygi committed
234 235 236 237 238 239 240
  localsize_.resize(nprow_);
  nrod_loc_.resize(nprow_);
  rod_h_.resize(nprow_);
  rod_k_.resize(nprow_);
  rod_lmin_.resize(nprow_);
  rod_size_.resize(nprow_);
  rod_first_.resize(nprow_);
241

Francois Gygi committed
242 243 244 245 246
  // resize with zero cutoff to initialize empty Basis
  resize(cell_,refcell_,0.0);
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
247
Basis::~Basis(void) {}
Francois Gygi committed
248 249

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
250
bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
Francois Gygi committed
251 252 253 254 255
  double ecut)
{
  assert(ecut>=0.0);
  assert(cell.volume() >= 0.0);
  assert(refcell.volume() >= 0.0);
256

257 258 259 260 261 262 263
  if ( ecut == ecut_ && refcell == refcell_ && refcell_.volume() != 0.0 )
  {
    cell_ = cell;
    // only the cell changes, ecut and the refcell remain unchanged
    update_g();
    return true;
  }
264

Francois Gygi committed
265 266 267
  ecut_ = ecut;
  cell_ = cell;
  refcell_ = refcell;
268

Francois Gygi committed
269 270 271 272 273 274 275 276
  if ( ecut == 0.0 || cell.volume() == 0.0)
  {
    idxmax_[0] = 0;
    idxmax_[1] = 0;
    idxmax_[2] = 0;
    idxmin_[0] = 0;
    idxmin_[1] = 0;
    idxmin_[2] = 0;
277

Francois Gygi committed
278 279 280 281 282 283 284 285 286 287 288 289 290
    size_ = 0;
    nrods_ = 0;
    for ( int ipe = 0; ipe < nprow_; ipe++ )
    {
      localsize_[ipe] = 0;
      nrod_loc_[ipe] = 0;
    }
    maxlocalsize_ = minlocalsize_ = 0;
    np_[0] = np_[1] = np_[2] = 0;
    idx_.resize(3*localsize_[myrow_]);
    g_.resize(localsize_[myrow_]);
    kpg_.resize(localsize_[myrow_]);
    gi_.resize(localsize_[myrow_]);
Francois Gygi committed
291
    kpgi_.resize(localsize_[myrow_]);
Francois Gygi committed
292 293 294
    g2_.resize(localsize_[myrow_]);
    kpg2_.resize(localsize_[myrow_]);
    g2i_.resize(localsize_[myrow_]);
Francois Gygi committed
295
    kpg2i_.resize(localsize_[myrow_]);
Francois Gygi committed
296
    gx_.resize(3*localsize_[myrow_]);
Francois Gygi committed
297
    kpgx_.resize(3*localsize_[myrow_]);
Francois Gygi committed
298 299 300
    isort_loc.resize(localsize_[myrow_]);
    return true;
  }
301

Francois Gygi committed
302 303 304 305 306 307
  const double two_ecut = 2.0 * ecut;
  const double twopi = 2.0 * M_PI;

  const double kpx = kpoint_.x;
  const double kpy = kpoint_.y;
  const double kpz = kpoint_.z;
308

Francois Gygi committed
309 310 311 312 313 314 315 316 317 318 319 320
  UnitCell defcell;
  // defcell: cell used to define which vectors are contained in the Basis
  // if refcell is defined, defcell = refcell
  // otherwise, defcell = cell
  if ( norm(refcell.b(0)) + norm(refcell.b(1)) + norm(refcell.b(2)) == 0.0 )
  {
    defcell = cell;
  }
  else
  {
    defcell = refcell;
  }
321

Francois Gygi committed
322 323 324
  const D3vector b0 = defcell.b(0);
  const D3vector b1 = defcell.b(1);
  const D3vector b2 = defcell.b(2);
325

Francois Gygi committed
326 327
  const double normb2 = norm(b2);
  const double b2inv2 = 1.0 / normb2;
328

Francois Gygi committed
329
  const D3vector kp = kpx*b0 + kpy*b1 + kpz*b2;
330

Francois Gygi committed
331
  if ( !cell.in_bz(kp) )
332 333 334 335 336
  {
    if ( ctxt_.onpe0() )
      cout << " Basis::resize: warning: " << kpoint_
           << " out of the BZ: " << kp << endl;
  }
Francois Gygi committed
337 338 339

  const double fac = sqrt(two_ecut) / twopi;

340
  // define safe enclosing domain for any k-point value in the BZ
341
  const int hmax = (int) ( 1.5 + fac * ( length(defcell.a(0) ) ) );
342
  const int hmin = - hmax;
343

344
  const int kmax = (int) ( 1.5 + fac * ( length(defcell.a(1) ) ) );
345
  const int kmin = - kmax;
346

347
  const int lmax = (int) ( 1.5 + fac * ( length(defcell.a(2) ) ) );
348
  const int lmin = - lmax;
349

Francois Gygi committed
350
  multiset<Rod> rodset;
351

Francois Gygi committed
352
  // build rod set
353

354 355 356 357 358 359
  int hmax_used = hmin;
  int hmin_used = hmax;
  int kmax_used = kmin;
  int kmin_used = kmax;
  int lmax_used = lmin;
  int lmin_used = lmax;
Francois Gygi committed
360 361 362 363 364 365

  if ( real_ )
  {
    // build basis at kpoint (0,0,0)

    // rod(0,0,0)
366 367 368 369
    // length of rod(0,0,0) is lend+1
    int lend = (int) ( sqrt(two_ecut * b2inv2) );
    size_ = lend + 1;
    rodset.insert(Rod(0,0,0,lend+1));
Francois Gygi committed
370 371
    nrods_ = 1;

372
    hmax_used = 0;
Francois Gygi committed
373
    hmin_used = 0;
374 375 376 377
    kmin_used = 0;
    kmax_used = 0;
    lmin_used = 0;
    lmax_used = lend;
378

Francois Gygi committed
379
    // rods (0,k,l)
380
    for ( int k = 1; k <= kmax; k++ )
Francois Gygi committed
381
    {
382 383
      int lstart=lmax,lend=lmin;
      bool found = false;
384
      for ( int l = lmin; l <= lmax; l++ )
Francois Gygi committed
385
      {
386 387
        const double two_e = norm(k*b1+l*b2);
        if ( two_e < two_ecut )
Francois Gygi committed
388
        {
389 390 391
          lstart = min(l,lstart);
          lend = max(l,lend);
          found = true;
Francois Gygi committed
392 393
        }
      }
394 395 396 397 398 399 400 401 402 403 404 405
      if ( found )
      {
        // non-zero intersection was found
        const int rodsize = lend - lstart + 1;
        size_ += rodsize;
        rodset.insert(Rod(0,k,lstart,rodsize));
        nrods_++;
        kmax_used = max(k,kmax_used);
        kmin_used = min(k,kmin_used);
        lmax_used = max(lend,lmax_used);
        lmin_used = min(lstart,lmin_used);
      }
Francois Gygi committed
406 407
    }
    // rods (h,k,l) h>0
408
    for ( int h = 1; h <= hmax; h++ )
Francois Gygi committed
409
    {
410
      for ( int k = kmin; k <= kmax; k++ )
Francois Gygi committed
411
      {
412 413
        int lstart=lmax,lend=lmin;
        bool found = false;
414
        for ( int l = lmin; l <= lmax; l++ )
Francois Gygi committed
415
        {
416 417
          const double two_e = norm(h*b0+k*b1+l*b2);
          if ( two_e < two_ecut )
Francois Gygi committed
418
          {
419 420 421
            lstart = min(l,lstart);
            lend = max(l,lend);
            found = true;
Francois Gygi committed
422 423
          }
        }
424 425 426 427 428 429 430 431 432 433 434 435 436 437
        if ( found )
        {
          // non-zero intersection was found
          const int rodsize = lend - lstart + 1;
          size_ += rodsize;
          rodset.insert(Rod(h,k,lstart,rodsize));
          nrods_++;
          hmax_used = max(h,hmax_used);
          hmin_used = min(h,hmin_used);
          kmax_used = max(k,kmax_used);
          kmin_used = min(k,kmin_used);
          lmax_used = max(lend,lmax_used);
          lmin_used = min(lstart,lmin_used);
        }
Francois Gygi committed
438 439 440 441 442 443 444 445 446 447
      }
    }
  }
  else
  {
    // build Basis for k != (0,0,0)

    size_ = 0;
    nrods_ = 0;
    // rods (h,k,l)
448
    for ( int h = hmin; h <= hmax; h++ )
Francois Gygi committed
449
    {
450
      for ( int k = kmin; k <= kmax; k++ )
Francois Gygi committed
451
      {
452 453
        int lstart=lmax,lend=lmin;
        bool found = false;
454
        for ( int l = lmin; l <= lmax; l++ )
Francois Gygi committed
455
        {
456 457
          const double two_e = norm((kpx+h)*b0 + (kpy+k)*b1 + (kpz+l)*b2);
          if ( two_e < two_ecut )
Francois Gygi committed
458
          {
459 460 461
            lstart = min(l,lstart);
            lend = max(l,lend);
            found = true;
Francois Gygi committed
462 463
          }
        }
464 465 466 467 468 469 470 471 472 473 474 475 476 477
        if ( found )
        {
          // non-zero intersection was found
          const int rodsize = lend - lstart + 1;
          size_ += rodsize;
          rodset.insert(Rod(h,k,lstart,rodsize));
          nrods_++;
          hmax_used = max(h,hmax_used);
          hmin_used = min(h,hmin_used);
          kmax_used = max(k,kmax_used);
          kmin_used = min(k,kmin_used);
          lmax_used = max(lend,lmax_used);
          lmin_used = min(lstart,lmin_used);
        }
Francois Gygi committed
478 479 480
      }
    }
  }
481

482 483 484 485 486 487 488 489
#if DEBUG
  cout << " hmin/hmax: " << hmin << " / " << hmax << endl;
  cout << " kmin/kmax: " << kmin << " / " << kmax << endl;
  cout << " lmin/lmax: " << lmin << " / " << lmax << endl;
  cout << " hmin/hmax used: " << hmin_used << " / " << hmax_used << endl;
  cout << " kmin/kmax used: " << kmin_used << " / " << kmax_used << endl;
  cout << " lmin/lmax used: " << lmin_used << " / " << lmax_used << endl;
#endif
Francois Gygi committed
490 491 492

  idxmax_[0] = hmax_used;
  idxmin_[0] = hmin_used;
493

Francois Gygi committed
494 495
  idxmax_[1] = kmax_used;
  idxmin_[1] = kmin_used;
496

Francois Gygi committed
497 498
  idxmax_[2] = lmax_used;
  idxmin_[2] = lmin_used;
499

500 501 502
  assert(hmax_used - hmin_used + 1 <= 2 * hmax);
  assert(kmax_used - kmin_used + 1 <= 2 * kmax);
  assert(lmax_used - lmin_used + 1 <= 2 * lmax);
503

Francois Gygi committed
504
  // compute good FFT sizes
505 506
  // use values independent of the kpoint
  int n;
507
  n = 2 * hmax;
508 509
  while ( !factorizable(n) ) n+=2;
  np_[0] = n;
510
  n = 2 * kmax;
511 512
  while ( !factorizable(n) ) n+=2;
  np_[1] = n;
513
  n = 2 * lmax;
514 515
  while ( !factorizable(n) ) n+=2;
  np_[2] = n;
Francois Gygi committed
516 517 518 519

  // Distribute the basis on nprow_ processors

  // build a min-heap of Nodes
520

Francois Gygi committed
521
  vector<Node*> nodes(nprow_);
522

Francois Gygi committed
523 524 525 526 527 528 529 530 531 532 533 534
  for ( int ipe = 0; ipe < nprow_; ipe++ )
  {
    nodes[ipe] = new Node(ipe);
    localsize_[ipe] = 0;
    nrod_loc_[ipe] = 0;
    rod_h_[ipe].resize(0);
    rod_k_[ipe].resize(0);
    rod_lmin_[ipe].resize(0);
    rod_size_[ipe].resize(0);
  }

  // nodes contains a valid min-heap of zero-size Nodes
535

Francois Gygi committed
536 537 538 539 540 541 542 543
  // insert rods into the min-heap
  // keep track of where rod(0,0,0) goes
  int pe_rod0 = -1, rank_rod0 = -1;
  multiset<Rod>::iterator p = rodset.begin();
  while ( p != rodset.end() )
  {
    // pop smallest element
    pop_heap(nodes.begin(), nodes.end(), ptr_greater<Node>());
544

Francois Gygi committed
545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
    // add rod size to smaller element
    nodes[nprow_-1]->addrod(*p);
    int ipe = nodes[nprow_-1]->id();

    // update info about rod on process ipe
    nrod_loc_[ipe]++;
    rod_h_[ipe].push_back(p->h());
    rod_k_[ipe].push_back(p->k());
    rod_lmin_[ipe].push_back(p->lmin());
    rod_size_[ipe].push_back(p->size());
    localsize_[ipe] += p->size();
    if ( p->h() == 0 && p->k() == 0 )
    {
      pe_rod0 = ipe;
      rank_rod0 = nodes[nprow_-1]->nrods()-1;
    }

    // push modified element back in the heap
    push_heap(nodes.begin(), nodes.end(), ptr_greater<Node>());
564

Francois Gygi committed
565 566
    p++;
  }
567 568

  maxlocalsize_ = (*max_element(nodes.begin(), nodes.end(),
Francois Gygi committed
569
    ptr_less<Node>()))->size();
570
  minlocalsize_ = (*min_element(nodes.begin(), nodes.end(),
Francois Gygi committed
571
    ptr_less<Node>()))->size();
572

Francois Gygi committed
573 574 575 576
  for ( int ipe = 0; ipe < nprow_; ipe++ )
  {
    delete nodes[ipe];
  }
577

Francois Gygi committed
578 579 580 581 582 583 584
  // swap node pe_rod0 with node 0 in order to have rod(0,0,0) on node 0
  swap(nrod_loc_[0], nrod_loc_[pe_rod0]);
  rod_h_[pe_rod0].swap(rod_h_[0]);
  rod_k_[pe_rod0].swap(rod_k_[0]);
  rod_lmin_[pe_rod0].swap(rod_lmin_[0]);
  rod_size_[pe_rod0].swap(rod_size_[0]);
  swap(localsize_[0], localsize_[pe_rod0]);
585
  //Node *tmpnodeptr = nodes[0]; nodes[0] = nodes[pe_rod0];
Francois Gygi committed
586
  //  nodes[pe_rod0]=tmpnodeptr;
587

Francois Gygi committed
588 589 590 591 592
  // reorder rods on node 0 so that rod(0,0,0) comes first
  swap(rod_h_[0][rank_rod0], rod_h_[0][0]);
  swap(rod_k_[0][rank_rod0], rod_k_[0][0]);
  swap(rod_lmin_[0][rank_rod0], rod_lmin_[0][0]);
  swap(rod_size_[0][rank_rod0], rod_size_[0][0]);
593

Francois Gygi committed
594 595 596 597 598 599 600 601 602 603 604
  // compute position of first element of rod (ipe,irod)
  for ( int ipe = 0; ipe < nprow_; ipe++ )
  {
    rod_first_[ipe].resize(nrod_loc_[ipe]);
    if ( nrod_loc_[ipe] > 0 )
      rod_first_[ipe][0] = 0;
    for ( int irod = 1; irod < nrod_loc_[ipe]; irod++ )
    {
      rod_first_[ipe][irod] = rod_first_[ipe][irod-1] + rod_size_[ipe][irod-1];
    }
  }
605

Francois Gygi committed
606
  // local arrays idx, g, gi, g2i, g2, gx
Francois Gygi committed
607 608 609 610 611 612 613 614 615
  idx_.resize(3*localsize_[myrow_]);
  int i = 0;
  for ( int irod = 0; irod < nrod_loc_[myrow_]; irod++ )
  {
    for ( int l = 0; l < rod_size_[myrow_][irod]; l++ )
    {
      idx_[3*i]   = rod_h_[myrow_][irod];
      idx_[3*i+1] = rod_k_[myrow_][irod];
      idx_[3*i+2] = rod_lmin_[myrow_][irod] + l;
616

Francois Gygi committed
617 618 619
      i++;
    }
  }
620

Francois Gygi committed
621 622 623
  g_.resize(localsize_[myrow_]);
  kpg_.resize(localsize_[myrow_]);
  gi_.resize(localsize_[myrow_]);
Francois Gygi committed
624
  kpgi_.resize(localsize_[myrow_]);
Francois Gygi committed
625 626 627
  g2_.resize(localsize_[myrow_]);
  kpg2_.resize(localsize_[myrow_]);
  g2i_.resize(localsize_[myrow_]);
Francois Gygi committed
628
  kpg2i_.resize(localsize_[myrow_]);
Francois Gygi committed
629
  gx_.resize(3*localsize_[myrow_]);
Francois Gygi committed
630
  kpgx_.resize(3*localsize_[myrow_]);
Francois Gygi committed
631
  isort_loc.resize(localsize_[myrow_]);
632

Francois Gygi committed
633
  update_g();
634

Francois Gygi committed
635
  // basis set construction is complete
636

Francois Gygi committed
637 638 639 640
  return true;
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
641
void Basis::update_g(void)
Francois Gygi committed
642
{
Francois Gygi committed
643
  // compute the values of g, kpg, gi, g2i, g2, kpg2, gx
Francois Gygi committed
644
  // N.B. use the values of cell (not defcell)
645

Francois Gygi committed
646 647 648 649 650 651
  const int locsize = localsize_[myrow_];
  for ( int i = 0; i < locsize; i++ )
  {
    D3vector gt = idx_[3*i+0] * cell_.b(0) +
                  idx_[3*i+1] * cell_.b(1) +
                  idx_[3*i+2] * cell_.b(2);
652

Francois Gygi committed
653 654 655
    D3vector kpgt = (kpoint_.x + idx_[3*i+0]) * cell_.b(0) +
                    (kpoint_.y + idx_[3*i+1]) * cell_.b(1) +
                    (kpoint_.z + idx_[3*i+2]) * cell_.b(2);
656

Francois Gygi committed
657 658 659
    gx_[i] = gt.x;
    gx_[locsize+i] = gt.y;
    gx_[locsize+locsize+i] = gt.z;
Francois Gygi committed
660 661 662
    kpgx_[i] = kpgt.x;
    kpgx_[locsize+i] = kpgt.y;
    kpgx_[locsize+locsize+i] = kpgt.z;
Francois Gygi committed
663 664 665

    g2_[i] = norm(gt);
    g_[i] = sqrt( g2_[i] );
666

Francois Gygi committed
667 668
    kpg2_[i] = norm(kpgt);
    kpg_[i] = sqrt( kpg2_[i] );
669

Francois Gygi committed
670
    gi_[i] = g_[i] > 0.0 ? 1.0 / g_[i] : 0.0;
Francois Gygi committed
671
    kpgi_[i] = kpg_[i] > 0.0 ? 1.0 / kpg_[i] : 0.0;
Francois Gygi committed
672
    g2i_[i] = gi_[i] * gi_[i];
Francois Gygi committed
673
    kpg2i_[i] = kpgi_[i] * kpgi_[i];
Francois Gygi committed
674 675
    isort_loc[i] = i;
  }
676

Francois Gygi committed
677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696
  VectorLess<double> g2_less(g2_);
  sort(isort_loc.begin(), isort_loc.end(), g2_less);
#if DEBUG
  for ( int i = 0; i < locsize; i++ )
  {
    cout << ctxt_.mype() << " sorted " << i << " " << g2_[isort_loc[i]] << endl;
  }
#endif
}

////////////////////////////////////////////////////////////////////////////////
void Basis::print(ostream& os)
{
  os << context().mype() << ": ";
  os << " Basis.kpoint():   " << kpoint() << endl;
  os << context().mype() << ": ";
  os << " Basis.kpoint():   " << kpoint().x << " * b0 + "
                              << kpoint().y << " * b1 + "
                              << kpoint().z << " * b2" << endl;
  os << context().mype() << ": ";
697
  os << " Basis.kpoint():   " << kpoint().x * cell().b(0) +
Francois Gygi committed
698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723
                                 kpoint().y * cell().b(1) +
                                 kpoint().z * cell().b(2) << endl;
  os << context().mype() << ": ";
  os << " Basis.cell():     " << endl << cell() << endl;
  os << context().mype() << ": ";
  os << " Basis.ref cell(): " << endl << refcell() << endl;
  os << context().mype() << ": ";
  os << " Basis.ecut():     " << ecut() << endl;
  os << context().mype() << ": ";
  os << " Basis.np(0,1,2):  " << np(0) << " "
       << np(1) << " " << np(2) << endl;
  os << context().mype() << ": ";
  os << " Basis.idxmin:       " << idxmin(0) << " "
       << idxmin(1) << " " << idxmin(2) << endl;
  os << context().mype() << ": ";
  os << " Basis.idxmax:       " << idxmax(0) << " "
       << idxmax(1) << " " << idxmax(2) << endl;
  os << context().mype() << ": ";
  os << " Basis.size():     " << size() << endl;
  os << context().mype() << ": ";
  os << " Basis.localsize(): " << localsize() << endl;
  os << context().mype() << ": ";
  os << " Basis.nrods():    " << nrods() << endl;
  os << context().mype() << ": ";
  os << " Basis.real():     " << real() << endl;
  os << context().mype() << ": ";
Francois Gygi committed
724
  os << " Basis total mem size (MB): " << memsize() / 1048576 << endl;
Francois Gygi committed
725
  os << context().mype() << ": ";
Francois Gygi committed
726
  os << " Basis local mem size (MB): " << localmemsize() / 1048576 << endl;
727

Francois Gygi committed
728
  os << context().mype() << ": ";
729
  os << "   ig      i   j   k        gx      gy      gz       |k+g|^2"
Francois Gygi committed
730 731
     << endl;
  os << context().mype() << ": ";
732
  os << "   --      -   -   -        --      --      --       -------"
Francois Gygi committed
733 734 735 736
     << endl;
  for ( int i = 0; i < localsize(); i++ )
  {
    os << context().mype() << ": ";
737 738 739
    os << setw(5) << i << "   "
       << setw(4) << idx(3*i)
       << setw(4) << idx(3*i+1)
Francois Gygi committed
740 741
       << setw(4) << idx(3*i+2)
       << "    "
Francois Gygi committed
742 743 744
       << setw(8) << setprecision(4) << gx(i)
       << setw(8) << setprecision(4) << gx(i+localsize())
       << setw(8) << setprecision(4) << gx(i+2*localsize())
Francois Gygi committed
745 746 747 748 749 750 751 752 753 754 755
       << setw(12) << setprecision(4) << 0.5 * kpg2(i)
       << endl;
  }
}

////////////////////////////////////////////////////////////////////////////////
ostream& operator<<(ostream& os, Basis& b)
{
  b.print(os);
  return os;
}