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 68
  // cout << " WavefunctionHandler::startElement " << StrX(qname) << endl;
  string locname(XMLString::transcode(localname));
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 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
  if ( locname == "wavefunction" || locname == "wavefunction_velocity" )
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      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 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
  }
  else if ( locname == "domain")
  {
    D3vector a,b,c;
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
  }
  else if ( locname == "reference_domain")
  {
    D3vector a,b,c;
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
  }
  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++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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 198 199 200 201 202
  }
  else if ( locname == "grid")
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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 243 244 245 246 247 248
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
    wftmp.resize((ft->np012loc()));
  }
  else if ( locname == "grid_function")
  {
    unsigned int len = attributes.getLength();
    for (unsigned int index = 0; index < len; index++)
    {
      string attrname(XMLString::transcode(attributes.getLocalName(index)));
      string attrval(XMLString::transcode(attributes.getValue(index)));
      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();
Francois Gygi committed
341 342 343 344 345
  string locname(XMLString::transcode(localname));
  //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 498 499 500 501
    const StructureHandler* const last)
{
  string locname(XMLString::transcode(localname));
  //cout << " WavefunctionHandler::endSubHandler " << locname << endl;
  delete last;
}