SampleReader.C 19.6 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// SampleReader.C:
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: SampleReader.C,v 1.28 2008-06-18 03:40:53 fgygi Exp $
Francois Gygi committed
7 8 9 10 11


#include "Sample.h"
#include "SampleReader.h"
#include "SpeciesReader.h"
Francois Gygi committed
12
#include "Species.h"
Francois Gygi committed
13 14 15 16 17 18 19 20 21 22 23
#include "Basis.h"
#include "FourierTransform.h"
#include "SlaterDet.h"

#include "XMLGFPreprocessor.h"

#include "Timer.h"

#include <cassert>
#include <string>
#include <iostream>
Francois Gygi committed
24
#include <fstream>
Francois Gygi committed
25 26
#include <sys/stat.h>

Francois Gygi committed
27 28 29
#if USE_XERCES
#include "SampleHandler.h"
#include "StructuredDocumentHandler.h"
Francois Gygi committed
30 31 32 33 34 35 36 37
#include <xercesc/util/XMLUniDefs.hpp>
#include <xercesc/sax2/Attributes.hpp>
#include <xercesc/util/PlatformUtils.hpp>
#include <xercesc/sax2/SAX2XMLReader.hpp>
#include <xercesc/framework/LocalFileInputSource.hpp>
#include <xercesc/sax2/XMLReaderFactory.hpp>
#include <xercesc/framework/MemBufInputSource.hpp>
using namespace xercesc;
Francois Gygi committed
38
#endif
39
using namespace std;
Francois Gygi committed
40 41 42 43 44 45 46

////////////////////////////////////////////////////////////////////////////////
SampleReader::SampleReader(const Context& ctxt) : ctxt_(ctxt) {}

////////////////////////////////////////////////////////////////////////////////
void SampleReader::readSample (Sample& s, const string uri, bool serial)
{
47 48
  Timer tm;
  tm.start();
Francois Gygi committed
49
#if USE_XERCES
Francois Gygi committed
50
  const char* encodingName = "UTF-8";
51 52
  SAX2XMLReader::ValSchemes valScheme = SAX2XMLReader::Val_Auto;
  //SAX2XMLReader::ValSchemes valScheme = SAX2XMLReader::Val_Always;
Francois Gygi committed
53 54 55 56 57 58 59
  //SAX2XMLReader::ValSchemes valScheme = SAX2XMLReader::Val_Never;
  bool expandNamespaces = false;
  bool doNamespaces = true;
  bool doSchema = true;
  bool schemaFullChecking = true;
  bool namespacePrefixes = false;
  SAX2XMLReader* parser = 0;
60

61 62
  Wavefunction* wfvtmp = new Wavefunction(s.wf);
  Wavefunction* current_wf = &s.wf;
Francois Gygi committed
63
  vector<vector<vector<double> > > dmat;
64
  int nx, ny, nz; // size of <grid> in wavefunction
Francois Gygi committed
65 66 67 68
  const int nspin = 1;
  const int current_ispin = 0;
  int current_ikp;
  dmat.resize(nspin);
69

Francois Gygi committed
70
  MemBufInputSource* memBufIS = 0;
71

Francois Gygi committed
72 73 74 75 76 77 78 79 80 81 82
  // XMLPlatformUtils initialization on task 0 only
  int ierr = 0;
  if ( ctxt_.onpe0() )
  {
    try
    {
      XMLPlatformUtils::Initialize();
    }

    catch (const XMLException& toCatch)
    {
83 84
      cout << "  Sample::readSample: Error during XML initialization :\n"
           << StrX(toCatch.getMessage()) << endl;
Francois Gygi committed
85 86 87 88 89 90 91 92 93 94 95
      ierr = 1;
    }
    ctxt_.ibcast_send(1,1,&ierr,1);
  }
  else
  {
    ctxt_.ibcast_recv(1,1,&ierr,1,0,0);
  }
  //cout << ctxt_.mype() <<": SampleReader: ierr=" << ierr << endl;
  if ( ierr > 0 )
    throw SampleReaderException("error in XMLPlatformUtils::Initialize");
96

Francois Gygi committed
97 98 99
  // Determine if uri is a local file
  struct stat statbuf;
  bool read_from_file = !stat(uri.c_str(),&statbuf);
100

Francois Gygi committed
101 102
  // check for serial override
  read_from_file &= !serial;
103

Francois Gygi committed
104 105 106 107 108 109 110 111 112
  string xmlcontent;
  DoubleMatrix gfdata(ctxt_);
  if ( read_from_file )
  {
    const char* const filename = uri.c_str();
    if ( ctxt_.onpe0() )
      cout << " SampleReader: reading from file "
           << filename << " size: "
           << statbuf.st_size << endl;
113

Francois Gygi committed
114
    XMLGFPreprocessor xmlgfp;
115

Francois Gygi committed
116
    xmlgfp.process(filename,gfdata,xmlcontent);
117

Francois Gygi committed
118 119
    if ( ctxt_.onpe0() )
    {
Francois Gygi committed
120 121 122 123 124
      cout << " xmlcontent.size(): " << xmlcontent.size()
           << endl;
#if DEBUG
      cout << ctxt_.mype() << ": xmlcontent: " << xmlcontent << endl;
#endif
Francois Gygi committed
125 126 127 128
      memBufIS = new MemBufInputSource
        ( (const XMLByte*) &xmlcontent[0], xmlcontent.size(), "buf_id", false);
    }
  }
129

Francois Gygi committed
130 131 132
  bool read_wf = false;
  bool read_wfv = false;
  // initialize wavefunction_velocity in case it appears in the sample file
133

Francois Gygi committed
134 135 136
  assert(sizeof(event_type)==sizeof(int));
  // remove default kpoint in wf
  s.wf.del_kpoint(D3vector(0.0,0.0,0.0));
137
  wfvtmp->del_kpoint(D3vector(0.0,0.0,0.0));
Francois Gygi committed
138

Francois Gygi committed
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
  if ( ctxt_.onpe0() )
  {
    parser = XMLReaderFactory::createXMLReader();
    if (valScheme == SAX2XMLReader::Val_Auto)
    {
        parser->setFeature(XMLUni::fgSAX2CoreValidation, true);
        parser->setFeature(XMLUni::fgXercesDynamic, true);
    }

    if (valScheme == SAX2XMLReader::Val_Never)
    {
        parser->setFeature(XMLUni::fgSAX2CoreValidation, false);
    }

    if (valScheme == SAX2XMLReader::Val_Always)
    {
        parser->setFeature(XMLUni::fgSAX2CoreValidation, true);
        parser->setFeature(XMLUni::fgXercesDynamic, false);
    }

    parser->setFeature(XMLUni::fgSAX2CoreNameSpaces, doNamespaces);
    parser->setFeature(XMLUni::fgXercesSchema, doSchema);
    parser->setFeature(XMLUni::fgXercesSchemaFullChecking, schemaFullChecking);
    parser->setFeature(XMLUni::fgSAX2CoreNameSpacePrefixes, namespacePrefixes);

    int errorCount = 0;
165 166
    SampleHandler* s_handler =
      new SampleHandler(s,gfdata,nx,ny,nz,dmat,*wfvtmp);
167

Francois Gygi committed
168 169 170 171 172 173 174 175 176 177 178 179
    try
    {
      StructuredDocumentHandler handler(s_handler);
      parser->setContentHandler(&handler);
      parser->setErrorHandler(&handler);

      cout << " Starting XML parsing" << endl;
      if ( read_from_file )
        parser->parse(*memBufIS);
      else
        parser->parse(uri.c_str());
      cout << " XML parsing done" << endl;
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
      errorCount = parser->getErrorCount();
    }

    catch (const XMLException& toCatch)
    {
        cout << "\nAn error occurred\n  Error: "
             << StrX(toCatch.getMessage())
             << "\n" << endl;
        //XMLPlatformUtils::Terminate();
        //delete parser;
        //throw;
    }

    catch (const SAXParseException& e)
    {
        cout << "\na SAXParseException occurred in file "
             << StrX(e.getSystemId())
             << ", line " << e.getLineNumber()
             << ", char " << e.getColumnNumber()
             << "\n  Message: " << StrX(e.getMessage()) << endl;
        //XMLPlatformUtils::Terminate();
        //delete parser;
        //throw;
    }
205

Francois Gygi committed
206
    read_wf = s_handler->read_wf;
Francois Gygi committed
207 208
    if ( read_wf )
      cout << " wavefunction was read" << endl;
Francois Gygi committed
209
    read_wfv = s_handler->read_wfv;
Francois Gygi committed
210 211
    if ( read_wfv )
      cout << " wavefunction velocity was read" << endl;
212

213 214 215
    cout << " SampleReader::readSample: grid nx,ny,nz="
         << nx << " " << ny << " " << nz << endl;

Francois Gygi committed
216 217 218
    delete s_handler;
    delete parser;
    XMLPlatformUtils::Terminate();
219

Francois Gygi committed
220
    // parsing of sample is complete, send end of sample tag to tasks > 0
Francois Gygi committed
221 222
    event_type event = end;
    ctxt_.ibcast_send(1,1,(int*)&event,1);
223

Francois Gygi committed
224 225 226 227 228 229
  } // onpe0
  else
  {
    // tasks > 0
    // listen for Sample events
    // cout << ctxt_.mype() << ": listening..." << endl;
230

Francois Gygi committed
231
    bool done = false;
Francois Gygi committed
232
    event_type event = invalid;
Francois Gygi committed
233 234
    while ( !done )
    {
Francois Gygi committed
235 236
      ctxt_.ibcast_recv(1,1,(int*)&event,1,0,0);
      if ( event == end )
Francois Gygi committed
237
        done = true;
238 239 240 241 242 243 244 245
      else if ( event == unit_cell )
      {
        // unit_cell
        double buf[9];
        ctxt_.dbcast_recv(9,1,buf,9,0,0);
        D3vector a(buf[0],buf[1],buf[2]);
        D3vector b(buf[3],buf[4],buf[5]);
        D3vector c(buf[6],buf[7],buf[8]);
Francois Gygi committed
246
        s.atoms.set_cell(a,b,c);
247
      }
Francois Gygi committed
248
      else if ( event == species )
Francois Gygi committed
249 250 251 252 253 254 255 256 257
      {
        // Species
        string curr_name;
        ctxt_.string_bcast(curr_name,0);
        //cout << ctxt_.mype() << ": receiving species " << curr_name << endl;
        Species* sp = new Species(ctxt_,curr_name);
        SpeciesReader sp_reader(ctxt_);
        sp_reader.bcastSpecies(*sp);
        s.atoms.addSpecies(sp,curr_name);
258
      }
Francois Gygi committed
259
      else if ( event == atom )
Francois Gygi committed
260 261 262 263 264 265 266
      {
        // Atom
        string curr_atom_name,curr_atom_species;
        D3vector curr_position, curr_velocity;
        ctxt_.string_bcast(curr_atom_name,0);
        ctxt_.string_bcast(curr_atom_species,0);
        double buf[3];
267
        ctxt_.dbcast_recv(3,1,buf,3,0,0);
Francois Gygi committed
268
        curr_position = D3vector(buf[0],buf[1],buf[2]);
269
        ctxt_.dbcast_recv(3,1,buf,3,0,0);
Francois Gygi committed
270 271 272 273 274 275
        curr_velocity = D3vector(buf[0],buf[1],buf[2]);
        //cout << ctxt_.mype() << ": receiving atom " << curr_atom_name << endl;
        Atom* a = new Atom(curr_atom_name, curr_atom_species,
                           curr_position, curr_velocity);
        s.atoms.addAtom(a);
      }
Francois Gygi committed
276
      else if ( event == wavefunction )
Francois Gygi committed
277
      {
Francois Gygi committed
278
        current_ikp = 0;
Francois Gygi committed
279 280 281 282 283 284
        // Wavefunction
        read_wf = true;
        int nel,nspin,nempty;
        ctxt_.ibcast_recv(1,1,&nel,1,0,0);
        ctxt_.ibcast_recv(1,1,&nspin,1,0,0);
        ctxt_.ibcast_recv(1,1,&nempty,1,0,0);
Francois Gygi committed
285
        //cout << ctxt_.mype() << ": SampleReader: receiving wf nel=" << nel
Francois Gygi committed
286 287 288 289
        //     << " nspin=" << nspin << " nempty=" << nempty << endl;
        s.wf.set_nel(nel);
        s.wf.set_nspin(nspin);
        s.wf.set_nempty(nempty);
Francois Gygi committed
290 291
        //cout << ctxt_.mype() << ": SampleReader: nel=" << s.wf.nel()
        //     << " nst=" << s.wf.nst() << endl;
292

Francois Gygi committed
293
        // domain
294
        double buf[9];
295
        ctxt_.dbcast_recv(9,1,buf,9,0,0);
Francois Gygi committed
296
        D3vector a(buf[0],buf[1],buf[2]);
297 298
        D3vector b(buf[3],buf[4],buf[5]);
        D3vector c(buf[6],buf[7],buf[8]);
Francois Gygi committed
299
        UnitCell uc(a,b,c);
300

Francois Gygi committed
301
        // grid
302 303 304 305 306 307
        int ibuf[3];
        ctxt_.ibcast_recv(3,1,ibuf,3,0,0);
        nx = ibuf[0];
        ny = ibuf[1];
        nz = ibuf[2];

Francois Gygi committed
308 309 310
        // receive only computed ecut
        double ecut;
        ctxt_.dbcast_recv(1,1,&ecut,1,0,0);
311

312
        // reference_domain
313
        ctxt_.dbcast_recv(9,1,buf,9,0,0);
314 315 316 317
        D3vector ar(buf[0],buf[1],buf[2]);
        D3vector br(buf[3],buf[4],buf[5]);
        D3vector cr(buf[6],buf[7],buf[8]);
        UnitCell ruc(ar,br,cr);
318

319
        s.wf.resize(uc,ruc,ecut);
320

Francois Gygi committed
321
        //cout << ctxt_.mype() << ": wf resized, ecut=" << ecut << endl;
322

Francois Gygi committed
323
      }
Francois Gygi committed
324
      else if ( event == wavefunction_velocity )
Francois Gygi committed
325
      {
326
        current_wf = wfvtmp; // set ptr used by the slater det section
Francois Gygi committed
327
        current_ikp = 0;
328 329
        //cout << ctxt_.mype()
        //     << ": SampleReader received wavefunction_velocity" << endl;
Francois Gygi committed
330 331 332 333 334 335 336 337
        read_wfv = true;
        // Wavefunction velocity
        int nel,nspin,nempty;
        ctxt_.ibcast_recv(1,1,&nel,1,0,0);
        ctxt_.ibcast_recv(1,1,&nspin,1,0,0);
        ctxt_.ibcast_recv(1,1,&nempty,1,0,0);
        //cout << ctxt_.mype() << ": receiving wf nel=" << nel
        //     << " nspin=" << nspin << " nempty=" << nempty << endl;
338 339 340
        wfvtmp->set_nel(nel);
        wfvtmp->set_nspin(nspin);
        wfvtmp->set_nempty(nempty);
341

Francois Gygi committed
342
        // domain
343
        double buf[9];
344
        ctxt_.dbcast_recv(9,1,buf,9,0,0);
Francois Gygi committed
345
        D3vector a(buf[0],buf[1],buf[2]);
346 347
        D3vector b(buf[3],buf[4],buf[5]);
        D3vector c(buf[6],buf[7],buf[8]);
Francois Gygi committed
348
        UnitCell uc(a,b,c);
349

Francois Gygi committed
350
        // grid
351 352 353 354 355 356
        int ibuf[3];
        ctxt_.ibcast_recv(3,1,ibuf,3,0,0);
        assert(nx == ibuf[0]);
        assert(ny == ibuf[1]);
        assert(nz == ibuf[2]);

Francois Gygi committed
357 358 359
        // receive only computed ecut
        double ecut;
        ctxt_.dbcast_recv(1,1,&ecut,1,0,0);
360

361
        // reference_domain
362
        ctxt_.dbcast_recv(9,1,buf,9,0,0);
363 364 365 366
        D3vector ar(buf[0],buf[1],buf[2]);
        D3vector br(buf[3],buf[4],buf[5]);
        D3vector cr(buf[6],buf[7],buf[8]);
        UnitCell ruc(ar,br,cr);
367

368
        current_wf->resize(uc,ruc,ecut);
369

Francois Gygi committed
370
        //cout << ctxt_.mype() << ": wfv resized, ecut=" << ecut << endl;
371

Francois Gygi committed
372 373 374 375 376 377
      }
      else if ( event == slater_determinant )
      {
        // process SlaterDet
        // receive kpoint and weight
        double buf[4];
378
        ctxt_.dbcast_recv(4,1,buf,4,0,0);
379
        current_wf->add_kpoint(D3vector(buf[0],buf[1],buf[2]),buf[3]);
Francois Gygi committed
380 381

        // receive density_matrix
382
        vector<double> dmat_tmp(current_wf->nst());
383
        s.wf.context().dbcast_recv(s.wf.nst(),1,&dmat_tmp[0],s.wf.nst(),0,0);
Francois Gygi committed
384 385 386
        dmat[current_ispin].push_back(dmat_tmp);

        if ( !read_from_file )
Francois Gygi committed
387
        {
Francois Gygi committed
388
          // receive grid_functions
389
          SlaterDet* sd = current_wf->sd(current_ispin,current_ikp);
Francois Gygi committed
390
          const Basis& basis = sd->basis();
391
          FourierTransform ft(basis,nx,ny,nz);
Francois Gygi committed
392 393 394 395 396 397
          const int wftmpr_size = basis.real() ? ft.np012loc() :
            2*ft.np012loc();
          valarray<double> wftmpr(wftmpr_size);
          vector<complex<double> > wftmp(ft.np012loc());
          int size = -1;
          for ( int nloc = 0; nloc < sd->nstloc(); nloc++ )
Francois Gygi committed
398
          {
Francois Gygi committed
399 400 401 402 403 404 405 406
            // read grid_function nloc fragment
            //cout << sd->context();
            //cout << sd->context().mype()
            //     << ": wf receiving nloc=" << nloc << endl;
            sd->context().irecv(1,1,&size,1,0,0);
            //cout << sd->context().mype()
            //     << ": received size=" << size << endl;
            assert(size==wftmpr_size);
407
            sd->context().drecv(size,1,&wftmpr[0],size,0,0);
Francois Gygi committed
408 409 410 411 412 413
            //cout << sd->context().mype()
            //     << ": grid_function nloc=" << nloc
            //     << "received" << endl;

            // copy to complex array
            if ( basis.real() )
Francois Gygi committed
414
            {
Francois Gygi committed
415
              for ( int i = 0; i < size; i++ )
Francois Gygi committed
416
              {
Francois Gygi committed
417
                wftmp[i] = wftmpr[i];
Francois Gygi committed
418 419
              }
            }
Francois Gygi committed
420 421 422 423 424 425 426 427 428 429
            else
            {
              memcpy((void*)&wftmp[0],(void*)&wftmpr[0],wftmpr_size*
                sizeof(double));
            }

            ComplexMatrix& c = sd->c();
            ft.forward(&wftmp[0],c.valptr(c.mloc()*nloc));
            //cout << sd->context().mype()
            //     << ": grid_function read for nloc=" << nloc << endl;
Francois Gygi committed
430 431 432
          }
        }
      }
Francois Gygi committed
433
    } // while !done receiving events from node 0
Francois Gygi committed
434 435 436 437 438
  } // if-else onpe0

  // This part is now executing on all tasks
  if ( read_from_file )
  {
439 440 441 442 443 444 445
    // dmat may contain two sets of density matrices if wfv was read
    // Only the first set of density matrices is used
    if ( read_wfv )
      assert(dmat[0].size() == 2 * s.wf.nkp() );
    else
      assert(dmat[0].size() == s.wf.nkp() );

446 447 448 449
#if DEBUG
    ofstream gffile("gf.dat");
    gffile << gfdata;
#endif
Francois Gygi committed
450 451 452
    if ( read_wf )
    {
      // transfer data from the gfdata matrix to the SlaterDets
Francois Gygi committed
453 454 455 456 457 458 459 460
#if DEBUG
      cout << ctxt_.mype() << ": mapping gfdata matrix..."
           << endl;
      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;
#endif
461

Francois Gygi committed
462 463 464 465 466 467 468
      for ( int ispin = 0; ispin < s.wf.nspin(); ispin++ )
      {
        for ( int ikp = 0; ikp < s.wf.nkp(); ikp++ )
        {
          SlaterDet* sd = s.wf.sd(ispin,ikp);
#if DEBUG
          cout << "SampleReader: ikp=" << ikp << " nst = " << sd->nst() << endl;
469
          cout << "SampleReader: ikp=" << ikp << " dmat size = "
Francois Gygi committed
470 471
               << dmat[ispin][ikp].size() << endl;
#endif
472

Francois Gygi committed
473 474 475
          // copy density matrix information
          sd->set_occ(dmat[ispin][ikp]);
          const Basis& basis = sd->basis();
476
          FourierTransform ft(basis,nx,ny,nz);
477 478 479 480
#if DEBUG
          cout << ctxt_.mype() << ": ft.np012loc()=" << ft.np012loc() << endl;
          cout << ctxt_.mype() << ": ft.context(): " << ft.context();
#endif
Francois Gygi committed
481 482 483 484

          ComplexMatrix& c = sd->c();
          // copy wf values
          // Determine the size of the temporary real matrix wftmpr
485
          int wftmpr_size, wftmpr_block_size;
486
          if ( basis.real() )
Francois Gygi committed
487 488
          {
            wftmpr_size = ft.np012();
489
            wftmpr_block_size = ft.np012loc(0);
Francois Gygi committed
490 491 492 493
          }
          else
          {
            wftmpr_size = 2*ft.np012();
494
            wftmpr_block_size = 2*ft.np012loc(0);
Francois Gygi committed
495
          }
496 497 498 499 500
#if DEBUG
          cout << ctxt_.mype() << ": wftmpr_size: " << wftmpr_size << endl;
          cout << ctxt_.mype() << ": wftmpr_block_size: "
               << wftmpr_block_size << endl;
#endif
Francois Gygi committed
501
          DoubleMatrix wftmpr(sd->context(),wftmpr_size,sd->nst(),
502
            wftmpr_block_size,c.nb());
Francois Gygi committed
503

Francois Gygi committed
504
          wftmpr.getsub(gfdata,wftmpr_size,sd->nst(),0,ikp*sd->nst());
505

506
#if DEBUG
Francois Gygi committed
507 508 509 510
          // 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
511
#endif
512

Francois Gygi committed
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535
          vector<complex<double> > wftmp(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];
            }
            else
            {
              // 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]);
            }
            ft.forward(&wftmp[0],c.valptr(c.mloc()*nloc));
          }
          // cout << " c matrix: ikp=" << ikp << endl;
          // cout << c;
        }
Francois Gygi committed
536 537
      }
    } // if read_wf
538

Francois Gygi committed
539 540 541
    if ( read_wfv )
    {
      // transfer wfv data from the gfdata matrix to the SlaterDets
542
      //cout << ctxt_.mype() << ": mapping gfdata matrix..."
Francois Gygi committed
543 544 545 546 547
      //     << endl;
      //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;
548

Francois Gygi committed
549
      for ( int ispin = 0; ispin < s.wf.nspin(); ispin++ )
Francois Gygi committed
550
      {
Francois Gygi committed
551 552
        for ( int ikp = 0; ikp < s.wf.nkp(); ikp++ )
        {
553
          SlaterDet* sd = wfvtmp->sd(ispin,ikp);
Francois Gygi committed
554
          const Basis& basis = sd->basis();
555
          FourierTransform ft(basis,nx,ny,nz);
Francois Gygi committed
556 557 558 559 560 561
          //cout << ctxt_.mype() << ": ft.np012loc()=" << ft.np012loc() << endl;
          //cout << ctxt_.mype() << ": ft.context(): " << ft.context();

          ComplexMatrix& c = sd->c();
          // copy wf values
          // Determine the size of the temporary real matrix wftmpr
562
          int wftmpr_size, wftmpr_block_size;
Francois Gygi committed
563 564 565 566
          if ( basis.real() && ( ft.np012() == gfdata.m() ) )
          {
            // all functions are real
            wftmpr_size = ft.np012();
567
            wftmpr_block_size = ft.np012loc(0);
Francois Gygi committed
568 569 570 571 572 573
          }
          else
          {
            // there is either 1) a mix of real and complex functions
            // or 2) only complex functions
            wftmpr_size = 2*ft.np012();
574
            wftmpr_block_size = 2*ft.np012loc(0);
Francois Gygi committed
575 576
          }
          DoubleMatrix wftmpr(sd->context(),wftmpr_size,sd->nst(),
577
            wftmpr_block_size,c.nb());
Francois Gygi committed
578

579 580
          wftmpr.getsub(gfdata,wftmpr_size,sd->nst(),0,
            (ikp+s.wf.nkp())*sd->nst());
Francois Gygi committed
581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602

          vector<complex<double> > wftmp(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];
            }
            else
            {
              // 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]);
            }
            ft.forward(&wftmp[0],c.valptr(c.mloc()*nloc));
          }
        }
Francois Gygi committed
603 604 605
      }
    }
  }
606

607 608
  // force consistency of unit cell
  // copy wavefunction domain on atomset unit_cell
Francois Gygi committed
609
  s.atoms.set_cell(s.wf.cell());
610

Francois Gygi committed
611
  // check if wavefunction_velocity element was read, if not, delete wfvtmp
612 613 614 615 616
  if ( s.wfv != 0 )
  {
    delete s.wfv;
    s.wfv = 0;
  }
Francois Gygi committed
617 618
  if ( read_wfv )
  {
619
    s.wfv = wfvtmp;
620
  }
Francois Gygi committed
621 622 623 624
#else
  // USE_XERCES was not defined
  if ( ctxt_.onpe0() )
  {
625
    cout << "  SampleReader: could not read (parser not defined)"
Francois Gygi committed
626 627 628
         << endl;
  }
#endif
629 630
  tm.stop();
  if ( ctxt_.onpe0() )
631
    cout << " SampleReader: read time: " << tm.real() << " s" << endl;
Francois Gygi committed
632
}