WavefunctionHandler.C 15.4 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
// WavefunctionHandler.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18

Francois Gygi committed
19 20 21 22
#include "WavefunctionHandler.h"
#include "Wavefunction.h"
#include "FourierTransform.h"
#include "Timer.h"
Francois Gygi committed
23
#include "SampleReader.h"
Francois Gygi committed
24 25 26 27 28 29 30 31

#include "StrX.h"
// XML transcoding for loading grid_functions
#include <xercesc/util/Base64.hpp>
#include <xercesc/util/XMLString.hpp>
using namespace xercesc;
#include <iostream>
#include <cassert>
32
#include <cstring> // memcpy
Francois Gygi committed
33 34 35 36
#include <sstream>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
37
WavefunctionHandler::WavefunctionHandler(Wavefunction& wf,
38 39
  DoubleMatrix& gfdata, int& current_gfdata_pos ) : wf_(wf), gfdata_(gfdata),
  current_gfdata_pos_(current_gfdata_pos) {}
Francois Gygi committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

////////////////////////////////////////////////////////////////////////////////
WavefunctionHandler::~WavefunctionHandler() {}

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::byteswap_double(size_t n, double* x)
{
  if (n==0) return;
  unsigned char* c = (unsigned char*) x;
  while ( n-- > 0 )
  {
    unsigned char tmp;
    tmp = c[7]; c[7] = c[0]; c[0] = tmp;
    tmp = c[6]; c[6] = c[1]; c[1] = tmp;
    tmp = c[5]; c[5] = c[2]; c[2] = tmp;
    tmp = c[4]; c[4] = c[3]; c[3] = tmp;
56

Francois Gygi committed
57 58
    c+=8;
  }
59
}
Francois Gygi committed
60 61 62 63 64

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::startElement(const XMLCh* const uri,
  const XMLCh* const localname, const XMLCh* const qname,
  const Attributes& attributes)
65
{
66
  bool onpe0 = wf_.context().onpe0();
Francois Gygi committed
67
  // cout << " WavefunctionHandler::startElement " << StrX(qname) << endl;
68
  string locname = StrX(localname).localForm();
69

Francois Gygi committed
70
  int nspin=1, nel=0, nempty=0;
71

Francois Gygi committed
72
  // consider only elements that are dealt with directly by WavefunctionHandler
73

Francois Gygi committed
74 75 76 77 78
  if ( locname == "wavefunction" || locname == "wavefunction_velocity" )
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
79
      string attrname = StrX(attributes.getLocalName(index)).localForm();
Francois Gygi committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
      if ( attrname == "ecut")
      {
        ecut = atof(StrX(attributes.getValue(index)).localForm());
      }
      else if ( attrname == "nspin")
      {
        nspin = atoi(StrX(attributes.getValue(index)).localForm());
      }
      else if ( attrname == "nempty" )
      {
        nempty = atoi(StrX(attributes.getValue(index)).localForm());
      }
      else if ( attrname == "nel" )
      {
        nel = atoi(StrX(attributes.getValue(index)).localForm());
      }
    }
97

98 99 100 101
    if ( onpe0 )
      cout << " WavefunctionHandler::startElement: " << locname
           << " nspin=" << nspin << " nel=" << nel << " nempty=" << nempty
           << endl;
102

Francois Gygi committed
103 104
    current_ispin = 0;
    current_ikp = 0;
105

Francois Gygi committed
106 107 108
    wf_.set_nel(nel);
    wf_.set_nspin(nspin);
    wf_.set_nempty(nempty);
109 110 111

    // remove default kpoint k=0
    wf_.del_kpoint(D3vector(0.0,0.0,0.0));
Francois Gygi committed
112 113 114 115 116 117 118
  }
  else if ( locname == "domain")
  {
    D3vector a,b,c;
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
119 120
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
Francois Gygi committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134
      istringstream stst(attrval);
      if ( attrname == "a")
      {
        stst >> a;
      }
      else if ( attrname == "b" )
      {
        stst >> b;
      }
      else if ( attrname == "c" )
      {
        stst >> c;
      }
    }
135

Francois Gygi committed
136 137 138
    //cout << " WavefunctionHandler::startElement: domain" << endl;
    uc.set(a,b,c);
    //cout << uc;
139 140 141 142 143 144 145
  }
  else if ( locname == "reference_domain")
  {
    D3vector a,b,c;
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
146 147
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
148 149 150 151 152 153 154 155 156 157 158 159 160 161
      istringstream stst(attrval);
      if ( attrname == "a")
      {
        stst >> a;
      }
      else if ( attrname == "b" )
      {
        stst >> b;
      }
      else if ( attrname == "c" )
      {
        stst >> c;
      }
    }
162

163 164 165
    //cout << " WavefunctionHandler::startElement: reference_domain" << endl;
    ruc.set(a,b,c);
    //cout << ruc;
Francois Gygi committed
166 167 168 169 170 171 172 173
  }
  else if ( locname == "density_matrix")
  {
    unsigned int len = attributes.getLength();
    int dmat_size = 0;
    string dmat_form;
    for (unsigned int index = 0; index < len; index++)
    {
174 175
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
Francois Gygi committed
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
      istringstream stst(attrval);
      if ( attrname == "form")
      {
        stst >> dmat_form;
      }
      else if ( attrname == "size" )
      {
        stst >> dmat_size;
      }
    }
    if ( dmat_form != "diagonal" )
    {
      cout << "WavefunctionHandler: density_matrix must be diagonal" << endl;
      wf_.context().abort(1);
    }
191
    dmat_.resize(dmat_size);
Francois Gygi committed
192 193 194 195 196 197
  }
  else if ( locname == "grid")
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
198 199
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
Francois Gygi committed
200 201 202
      istringstream stst(attrval);
      if ( attrname == "nx")
      {
203
        stst >> nx_;
Francois Gygi committed
204 205 206
      }
      else if ( attrname == "ny" )
      {
207
        stst >> ny_;
Francois Gygi committed
208 209 210
      }
      else if ( attrname == "nz" )
      {
211
        stst >> nz_;
Francois Gygi committed
212 213
      }
    }
214

Francois Gygi committed
215 216 217 218
    if ( ecut == 0.0 )
    {
      // ecut attribute was not specified. Infer from grid size
      // Ecut = max(ecut_x,ecut_y,ecut_z)
219

Francois Gygi committed
220 221 222
      // When importing grids with Dirichlet b.c. grid sizes can be odd
      // round nx,ny,nz to next even number to compute ecut
      // use nx+nx%2 instead of nx
223 224 225
      double g0_max = ((2*(nx_+nx_%2)-1.0)/2.0) * 2.0 * M_PI / length(uc.a(0));
      double g1_max = ((2*(ny_+ny_%2)-1.0)/2.0) * 2.0 * M_PI / length(uc.a(1));
      double g2_max = ((2*(nz_+nz_%2)-1.0)/2.0) * 2.0 * M_PI / length(uc.a(2));
Francois Gygi committed
226 227 228 229 230 231 232
      double ecut0 = 0.125 * g0_max * g0_max;
      double ecut1 = 0.125 * g1_max * g1_max;
      double ecut2 = 0.125 * g2_max * g2_max;

      ecut = max(max(ecut0,ecut1),ecut2);
      cout << " ecut=" << 2*ecut << " Ry" << endl;
    }
233

234
    wf_.resize(uc,ruc,ecut);
Francois Gygi committed
235 236 237
  }
  else if ( locname == "slater_determinant")
  {
238 239
    if ( onpe0 )
      cout << " WavefunctionHandler::startElement: slater_determinant" << endl;
Francois Gygi committed
240 241 242
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
243 244
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
Francois Gygi committed
245 246 247 248
      istringstream stst(attrval);
      if ( attrname == "kpoint")
      {
        stst >> current_kx >> current_ky >> current_kz;
249 250 251
        if ( onpe0 )
          cout << " kpoint=" << current_kx
               << " " << current_ky << " " << current_kz;
Francois Gygi committed
252 253 254 255
      }
      else if ( attrname == "weight" )
      {
        stst >> current_weight;
256 257
        if ( onpe0 )
          cout << " weight=" << current_weight;
Francois Gygi committed
258 259 260 261
      }
      else if ( attrname == "size" )
      {
        stst >> current_size;
262 263
        if ( onpe0 )
          cout << " size=" << current_size;
Francois Gygi committed
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279
      }
      else if ( attrname == "spin" )
      {
        std::string spin;
        stst >> spin;
        if (spin == "up" ) current_ispin=0;
        else if (spin == "down" )
        {
          int last_ispin=current_ispin;
          current_ispin=1;
          // reset kpoint index if necessary
          if ( last_ispin==0 && current_ispin==1 )
          {
            current_ikp=0;
          }
        }
280 281
        if ( onpe0 )
          cout << " read slater_determinant: spin=" << spin;
Francois Gygi committed
282 283
      }
    }
284 285 286
    if ( onpe0 )
      cout << endl;

Francois Gygi committed
287 288 289 290
    // add kpoint only if spin is up (if spin is down,
    // the kpoint should be defined already)
    if ( current_ispin == 0 )
      wf_.add_kpoint(D3vector(current_kx,current_ky,current_kz),current_weight);
291 292 293 294

    // resize sd(current_ispin,current_ikp)->nst() == current_size
    wf_.sd(current_ispin,current_ikp)->resize(wf_.cell(),
      wf_.refcell(), wf_.ecut(), current_size);
Francois Gygi committed
295

Francois Gygi committed
296
    const Basis& basis = wf_.sd(current_ispin,current_ikp)->basis();
297
    ft = new FourierTransform(basis,nx_,ny_,nz_);
Francois Gygi committed
298 299 300 301 302 303 304
    wftmp.resize((ft->np012loc()));
  }
  else if ( locname == "grid_function")
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
305 306
      string attrname = StrX(attributes.getLocalName(index)).localForm();
      string attrval = StrX(attributes.getValue(index)).localForm();
Francois Gygi committed
307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
      istringstream stst(attrval);
      if ( attrname == "nx")
      {
        stst >> current_gf_nx;
      }
      else if ( attrname == "ny" )
      {
        stst >> current_gf_ny;
      }
      else if ( attrname == "nz" )
      {
        stst >> current_gf_nz;
      }
      else if ( attrname == "encoding" )
      {
        stst >> current_gf_encoding;
      }
    }
325

Francois Gygi committed
326
    //cout << " WavefunctionHandler::startElement: grid_function"
327 328
    //     << " nx=" << current_gf_nx
    //     << " ny=" << current_gf_ny
Francois Gygi committed
329 330 331
    //     << " nz=" << current_gf_nz
    //     << "\n encoding=" << current_gf_encoding
    //     << endl;
332
  }
Francois Gygi committed
333 334 335 336 337 338
}

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::endElement(const XMLCh* const uri,
  const XMLCh* const localname, const XMLCh* const qname, string& content)
{
339 340
  const Context& ctxt = wf_.context();
  bool onpe0 = ctxt.onpe0();
341
  string locname = StrX(localname).localForm();
Francois Gygi committed
342 343 344 345
  //cout << " WavefunctionHandler::endElement " << locname << endl;
  if ( locname == "density_matrix")
  {
    istringstream stst(content);
346 347
    for ( int i = 0; i < dmat_.size(); i++ )
      stst >> dmat_[i];
Francois Gygi committed
348 349 350
  }
  else if ( locname == "grid_function")
  {
351 352 353 354 355
    // current implementation accepts only full grids as declared in
    // the wavefunction <grid> element
    assert(current_gf_nx==nx_ &&
           current_gf_ny==ny_ &&
           current_gf_nz==nz_ );
356 357 358 359 360 361 362 363 364 365
  }
  else if ( locname == "slater_determinant")
  {
    // data was read into the gfdata matrix
    // transfer data from the gfdata matrix to the SlaterDet object
    if ( onpe0 )
      cout << " WavefunctionHandler::endElement: slater_determinant" << endl;
    SlaterDet* sd = wf_.sd(current_ispin,current_ikp);

#if DEBUG
366
    cout << ctxt.mype() << ": mapping gfdata matrix..."
367
         << endl;
368 369 370 371
    cout << ctxt.mype() << ": gfdata: (" << gfdata_.m() << "x" << gfdata_.n()
         << ") (" << gfdata_.mb() << "x" << gfdata_.nb() << ") blocks" << endl;
    cout << ctxt.mype() << ": gfdata.mloc()=" << gfdata_.mloc()
         << " gfdata.nloc()=" << gfdata_.nloc() << endl;
372 373 374 375 376
#endif

    sd->set_occ(dmat_);
    const Basis& basis = sd->basis();
#if DEBUG
377 378
    cout << ctxt.mype() << ": ft->np012loc()=" << ft->np012loc() << endl;
    cout << ctxt.mype() << ": ft->context(): " << ft->context();
379
#endif
380

381 382 383 384 385
    ComplexMatrix& c = sd->c();
    // copy wf values
    // Determine the size of the temporary real matrix wftmpr
    int wftmpr_size, wftmpr_block_size;
    if ( basis.real() )
Francois Gygi committed
386
    {
387 388
      wftmpr_size = ft->np012();
      wftmpr_block_size = ft->np012loc(0);
Francois Gygi committed
389 390 391
    }
    else
    {
392 393 394 395
      wftmpr_size = 2*ft->np012();
      wftmpr_block_size = 2*ft->np012loc(0);
    }
#if DEBUG
396 397
    cout << ctxt.mype() << ": wftmpr_size: " << wftmpr_size << endl;
    cout << ctxt.mype() << ": wftmpr_block_size: "
398 399 400 401
         << wftmpr_block_size << endl;
#endif
    DoubleMatrix wftmpr(sd->context(),wftmpr_size,sd->nst(),
                        wftmpr_block_size,c.nb());
Francois Gygi committed
402

403 404 405 406 407 408 409
#if DEBUG
    // parameters of the getsub operation
    cout << "WavefunctionHandler: current_ikp=  " << current_ikp << endl;
    cout << "WavefunctionHandler: current_ispin=" << current_ispin << endl;
    cout << "WavefunctionHandler: sd->nst()=    " << sd->nst() << endl;
    cout << "WavefunctionHandler: wf_.nkp()=    " << wf_.nkp() << endl;
    cout << "WavefunctionHandler: current_gfdata_pos= "
410
     << current_gfdata_pos_ << endl;
411 412
#endif

413 414 415
    assert(current_gfdata_pos_ < gfdata_.n());
    wftmpr.getsub(gfdata_,wftmpr_size,sd->nst(),0,current_gfdata_pos_);
    current_gfdata_pos_ += sd->nst();
Francois Gygi committed
416

417 418 419 420 421
#if DEBUG
    // Check orthogonality by computing overlap matrix
    DoubleMatrix smat(sd->context(),sd->nst(),sd->nst(),c.nb(),c.nb());
    smat.syrk('l','t',1.0,wftmpr,0.0);
    cout << smat;
Francois Gygi committed
422
#endif
423 424 425 426 427 428 429 430 431 432 433

    wftmp.resize(ft->np012loc());
    for ( int nloc = 0; nloc < sd->nstloc(); nloc++ )
    {
      // copy column of wftmpr to complex array wftmp
      if ( wftmpr_size == ft->np012() )
      {
        // function is real and must be expanded
        double* p = wftmpr.valptr(wftmpr.mloc()*nloc);
        for ( int i = 0; i < ft->np012loc(); i++ )
          wftmp[i] = p[i];
Francois Gygi committed
434 435 436
      }
      else
      {
437 438 439 440
        // function is complex
        double* p = wftmpr.valptr(wftmpr.mloc()*nloc);
        for ( int i = 0; i < ft->np012loc(); i++ )
          wftmp[i] = complex<double>(p[2*i],p[2*i+1]);
Francois Gygi committed
441
      }
442 443
      ft->forward(&wftmp[0],c.valptr(c.mloc()*nloc));
    }
Francois Gygi committed
444

445 446 447 448 449 450 451 452 453
    current_ikp++;
    delete ft;
  }
  else if ( locname == "wavefunction" || locname == "wavefunction_velocity" )
  {
    // check that all SlaterDets of a given spin and have same nst
    for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
    {
      for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
Francois Gygi committed
454
      {
455
        if ( wf_.sd_[ispin][ikp]->nst() != wf_.sd_[ispin][0]->nst() )
Francois Gygi committed
456
        {
457 458
          cout << "nst differs for different kpoints in sample file" << endl;
          wf_.ctxt_.abort(1);
Francois Gygi committed
459 460
        }
      }
461
    }
Francois Gygi committed
462

463 464 465 466
    // set nst_[ispin] to match the values read
    for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
      wf_.nst_[ispin] = wf_.sd_[ispin][0]->nst();

467
#if 1
468 469 470 471 472 473 474 475 476
    // if nspin==2, adjust deltaspin_ to reflect the number of states
    // of each spin that were read
    if ( wf_.nspin() == 2 )
    {
      // assume that up spin is the majority spin
      assert(wf_.nst(0) >= wf_.nst(1));
      // The following line is correct for both
      // even and odd number of electrons
      wf_.deltaspin_ = (wf_.nst(0)-wf_.nst(1))/2;
Francois Gygi committed
477
    }
478
#endif
Francois Gygi committed
479 480 481 482 483
  }
}

////////////////////////////////////////////////////////////////////////////////
StructureHandler* WavefunctionHandler::startSubHandler(const XMLCh* const uri,
484
    const XMLCh* const localname, const XMLCh* const qname,
Francois Gygi committed
485 486 487 488 489 490 491 492 493 494
    const Attributes& attributes)
{
  // check if element qname can be processed by another StructureHandler
  // If it can, return a pointer to the StructureHandler, otherwise return 0
  // cout << " WavefunctionHandler::startSubHandler " << StrX(qname) << endl;
  return 0;
}

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::endSubHandler(const XMLCh* const uri,
495
    const XMLCh* const localname, const XMLCh* const qname,
Francois Gygi committed
496 497
    const StructureHandler* const last)
{
498
  string locname = StrX(localname).localForm();
Francois Gygi committed
499 500 501
  //cout << " WavefunctionHandler::endSubHandler " << locname << endl;
  delete last;
}