WavefunctionHandler.C 14.7 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// WavefunctionHandler.C
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: WavefunctionHandler.C,v 1.15 2008-01-13 23:04:46 fgygi Exp $
Francois Gygi committed
7 8

#if USE_XERCES
Francois Gygi committed
9 10 11 12 13

#include "WavefunctionHandler.h"
#include "Wavefunction.h"
#include "FourierTransform.h"
#include "Timer.h"
Francois Gygi committed
14
#include "SampleReader.h"
Francois Gygi committed
15 16 17 18 19 20 21 22 23 24 25 26

#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>
#include <sstream>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
27
WavefunctionHandler::WavefunctionHandler(Wavefunction& wf,
Francois Gygi committed
28 29 30
  DoubleMatrix& gfdata,
  vector<vector<vector<double> > > &dmat) :
  wf_(wf), gfdata_(gfdata), dmat_(dmat), ecut(0.0)
Francois Gygi committed
31
{
Francois Gygi committed
32 33
  // if the data is read from a file, gfdata has a finite size
  // since the grid functions were processed by XMLGFPreprocessor
Francois Gygi committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  // if the gfdata matrix has finite dimensions, set read_from_gfdata flag
  read_from_gfdata = ( gfdata.n() > 0 );
  current_igf = 0;
}

////////////////////////////////////////////////////////////////////////////////
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;
54

Francois Gygi committed
55 56
    c+=8;
  }
57
}
Francois Gygi committed
58 59 60 61 62

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::startElement(const XMLCh* const uri,
  const XMLCh* const localname, const XMLCh* const qname,
  const Attributes& attributes)
63
{
Francois Gygi committed
64 65
  // cout << " WavefunctionHandler::startElement " << StrX(qname) << endl;
  string locname(XMLString::transcode(localname));
66

Francois Gygi committed
67
  int nspin=1, nel=0, nempty=0;
Francois Gygi committed
68
  event_type event = invalid;
Francois Gygi committed
69
  // consider only elements that are dealt with directly by WavefunctionHandler
70

Francois Gygi committed
71 72
  if ( locname == "wavefunction" || locname == "wavefunction_velocity" )
  {
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
    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());
      }
    }
95

Francois Gygi committed
96 97 98
    cout << " WavefunctionHandler::startElement: " << locname
         << " nspin=" << nspin << " nel=" << nel << " nempty=" << nempty
         << endl;
99

Francois Gygi committed
100 101 102
    current_ispin = 0;
    current_ikp = 0;
    current_n = 0;
103

Francois Gygi committed
104 105
    // notify listening nodes
    if ( locname == "wavefunction" )
Francois Gygi committed
106
      event = wavefunction;
Francois Gygi committed
107
    else
Francois Gygi committed
108
      event = wavefunction_velocity;
109

Francois Gygi committed
110
    wf_.context().ibcast_send(1,1,(int*)&event,1);
Francois Gygi committed
111 112 113
    wf_.context().ibcast_send(1,1,&nel,1);
    wf_.context().ibcast_send(1,1,&nspin,1);
    wf_.context().ibcast_send(1,1,&nempty,1);
114

Francois Gygi committed
115 116 117
    wf_.set_nel(nel);
    wf_.set_nspin(nspin);
    wf_.set_nempty(nempty);
Francois Gygi committed
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
    assert(nspin==1);
  }
  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;
      }
    }
142

Francois Gygi committed
143 144 145
    //cout << " WavefunctionHandler::startElement: domain" << endl;
    uc.set(a,b,c);
    //cout << uc;
146

Francois Gygi committed
147
    // notify listening nodes
148 149 150 151
    double buf[9];
    buf[0] = uc.a(0).x; buf[1] = uc.a(0).y; buf[2] = uc.a(0).z;
    buf[3] = uc.a(1).x; buf[4] = uc.a(1).y; buf[5] = uc.a(1).z;
    buf[6] = uc.a(2).x; buf[7] = uc.a(2).y; buf[8] = uc.a(2).z;
152
    wf_.context().dbcast_send(9,1,buf,9);
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
  }
  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;
      }
    }
176

177 178 179
    //cout << " WavefunctionHandler::startElement: reference_domain" << endl;
    ruc.set(a,b,c);
    //cout << ruc;
180

Francois Gygi committed
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
  }
  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);
    }
Francois Gygi committed
206
    dmat_tmp.resize(dmat_size);
Francois Gygi committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
  }
  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")
      {
        stst >> nx;
      }
      else if ( attrname == "ny" )
      {
        stst >> ny;
      }
      else if ( attrname == "nz" )
      {
        stst >> nz;
      }
    }
229

Francois Gygi committed
230 231 232 233
    if ( ecut == 0.0 )
    {
      // ecut attribute was not specified. Infer from grid size
      // Ecut = max(ecut_x,ecut_y,ecut_z)
234

Francois Gygi committed
235 236 237 238 239 240 241 242 243 244 245 246 247
      // 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
      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));
      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;
    }
248

Francois Gygi committed
249 250
    // notify listening nodes of ecut
    wf_.context().dbcast_send(1,1,&ecut,1);
251

252 253 254 255 256 257 258
    // notify listening nodes of the reference_domain
    // note: the reference_domain is optional in the sample file
    // notify listening nodes
    double buf[9];
    buf[0] = ruc.a(0).x; buf[1] = ruc.a(0).y; buf[2] = ruc.a(0).z;
    buf[3] = ruc.a(1).x; buf[4] = ruc.a(1).y; buf[5] = ruc.a(1).z;
    buf[6] = ruc.a(2).x; buf[7] = ruc.a(2).y; buf[8] = ruc.a(2).z;
259
    wf_.context().dbcast_send(9,1,buf,9);
260

261
    wf_.resize(uc,ruc,ecut);
Francois Gygi committed
262 263 264
  }
  else if ( locname == "slater_determinant")
  {
Francois Gygi committed
265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
    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;
        cout << " read <slater_determinant> kpoint=" << current_kx
             << " " << current_ky << " " << current_kz;
      }
      else if ( attrname == "weight" )
      {
        stst >> current_weight;
        cout << " weight=" << current_weight;
      }
      else if ( attrname == "size" )
      {
        stst >> current_size;
        cout << " size=" << current_size << endl;
      }
    }
    // notify listening nodes: slater_determinant
    event = slater_determinant;
    wf_.context().ibcast_send(1,1,(int*)&event,1);

    // send kpoint and weight to listening nodes
    double buf[4];
    buf[0] = current_kx; buf[1] = current_ky; buf[2] = current_kz;
    buf[3] = current_weight;
296
    wf_.context().dbcast_send(4,1,buf,4);
Francois Gygi committed
297 298 299
    wf_.add_kpoint(D3vector(current_kx,current_ky,current_kz),current_weight);
    assert(current_size==wf_.sd(current_ispin,current_ikp)->nst());

Francois Gygi committed
300
    const Basis& basis = wf_.sd(current_ispin,current_ikp)->basis();
Francois Gygi committed
301

Francois Gygi committed
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
    // check the size of the basis generated
    //cout << " sd basis: np0,np1,np2 = " << basis.np(0)
    //     << " " << basis.np(1)
    //     << " " << basis.np(2)
    //     << endl;
    ft = new FourierTransform(basis,basis.np(0),basis.np(1),basis.np(2));
    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;
      }
    }
335

Francois Gygi committed
336
    //cout << " WavefunctionHandler::startElement: grid_function"
337 338
    //     << " nx=" << current_gf_nx
    //     << " ny=" << current_gf_ny
Francois Gygi committed
339 340 341
    //     << " nz=" << current_gf_nz
    //     << "\n encoding=" << current_gf_encoding
    //     << endl;
342
  }
Francois Gygi committed
343 344 345 346 347 348 349 350 351 352 353
}

////////////////////////////////////////////////////////////////////////////////
void WavefunctionHandler::endElement(const XMLCh* const uri,
  const XMLCh* const localname, const XMLCh* const qname, string& content)
{
  string locname(XMLString::transcode(localname));
  //cout << " WavefunctionHandler::endElement " << locname << endl;
  if ( locname == "density_matrix")
  {
    istringstream stst(content);
Francois Gygi committed
354 355
    for ( int i = 0; i < dmat_tmp.size(); i++ )
      stst >> dmat_tmp[i];
356

Francois Gygi committed
357
    // send dmat to listening nodes
Francois Gygi committed
358
    dmat_[current_ispin].push_back(dmat_tmp);
359
    wf_.context().dbcast_send(dmat_tmp.size(),1,&dmat_tmp[0],dmat_tmp.size());
Francois Gygi committed
360 361 362 363 364 365 366
  }
  else if ( locname == "grid_function")
  {
    // current implementation accepts only full grids
    assert(current_gf_nx==ft->np0());
    assert(current_gf_ny==ft->np1());
    assert(current_gf_nz==ft->np2());
367

Francois Gygi committed
368 369 370 371 372 373 374 375 376
    if ( read_from_gfdata )
    {
      // do nothing
      //cout << "WavefunctionHandler::endElement: current_igf=" << current_igf
      //     << endl;
      current_igf++;
    }
    else
    {
Francois Gygi committed
377 378 379 380 381 382 383 384
      const Basis& basis = wf_.sd(current_ispin,current_ikp)->basis();
      int wftmpr_size = current_gf_nx*current_gf_ny*current_gf_nz;
      // if basis is complex, double the size of wftmpr
      if ( !basis.real() )
        wftmpr_size *= 2;

      valarray<double> wftmpr(wftmpr_size);

Francois Gygi committed
385 386 387
      if ( current_gf_encoding == "text" )
      {
        istringstream stst(content);
Francois Gygi committed
388 389
        for ( int i = 0; i < wftmpr_size; i++ )
          stst >> wftmpr[i];
Francois Gygi committed
390 391 392
      }
      else if ( current_gf_encoding == "base64" )
      {
Francois Gygi committed
393
        // base64 encoding
Francois Gygi committed
394 395 396 397
        unsigned int length;
        XMLByte* b = Base64::decode((XMLByte*)content.c_str(), &length);
        assert(b!=0);
        // use data in b
Francois Gygi committed
398
        assert(length/sizeof(double)==wftmpr_size);
Francois Gygi committed
399
#if PLT_BIG_ENDIAN
400
        byteswap_double(wftmpr_size,(double*)b);
Francois Gygi committed
401
#endif
Francois Gygi committed
402 403 404 405 406 407 408 409 410 411 412 413
        memcpy(&wftmpr[0],b,wftmpr_size*sizeof(double));
        XMLString::release(&b);
      }
      else
      {
        cout << "WavefunctionHandler: unknown grid_function encoding" << endl;
        return;
      }

      // wftmpr now contains the grid function

      // send subgrids to listening nodes
Francois Gygi committed
414

Francois Gygi committed
415 416 417 418 419 420 421 422 423
      SlaterDet* sd = wf_.sd(current_ispin,current_ikp);
      assert(sd != 0);
      // pcol = process column destination
      int pcol = sd->c().pc(current_n);
      for ( int prow = 0; prow < sd->context().nprow(); prow++ )
      {
        int size = ft->np2_loc(prow) * ft->np0() * ft->np1();
        int istart = ft->np2_first(prow) * ft->np0() * ft->np1();
        if ( !basis.real() )
Francois Gygi committed
424
        {
Francois Gygi committed
425 426
          size *= 2;
          istart *= 2;
Francois Gygi committed
427
        }
Francois Gygi committed
428 429
        // send subgrid to node (prow,pcol)
        if ( !(prow==0 && pcol==0) )
Francois Gygi committed
430
        {
Francois Gygi committed
431
          sd->context().isend(1,1,&size,1,prow,pcol);
432
          sd->context().dsend(size,1,&wftmpr[istart],size,prow,pcol);
Francois Gygi committed
433 434
        }
      }
Francois Gygi committed
435 436 437 438

      // if destination column is pcol=0, copy to complex array on node 0
      // and process
      if ( pcol == 0 )
Francois Gygi committed
439
      {
Francois Gygi committed
440 441 442 443 444 445 446 447 448 449 450 451
        if ( basis.real() )
        {
          for ( int i = 0; i < ft->np012loc(); i++ )
            wftmp[i] = wftmpr[i];
        }
        else
        {
          for ( int i = 0; i < ft->np012loc(); i++ )
            wftmp[i] = complex<double>(wftmpr[2*i],wftmpr[2*i+1]);
        }
        ComplexMatrix& c = wf_.sd(current_ispin,current_ikp)->c();
        ft->forward(&wftmp[0],c.valptr(c.mloc()*current_n));
Francois Gygi committed
452 453 454 455 456 457 458
      }
    }
    //cout << " grid_function read n=" << current_n << endl;
    current_n++;
  }
  else if ( locname == "slater_determinant")
  {
Francois Gygi committed
459
    current_ikp++;
Francois Gygi committed
460 461 462 463 464 465
    delete ft;
  }
}

////////////////////////////////////////////////////////////////////////////////
StructureHandler* WavefunctionHandler::startSubHandler(const XMLCh* const uri,
466
    const XMLCh* const localname, const XMLCh* const qname,
Francois Gygi committed
467 468 469 470 471 472 473 474 475 476
    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,
477
    const XMLCh* const localname, const XMLCh* const qname,
Francois Gygi committed
478 479 480 481 482 483
    const StructureHandler* const last)
{
  string locname(XMLString::transcode(localname));
  //cout << " WavefunctionHandler::endSubHandler " << locname << endl;
  delete last;
}
484

Francois Gygi committed
485
#endif