Fixed compatibility problem with older version for reading samples with

smaller grid sizes.


git-svn-id: http://qboxcode.org/svn/qb/trunk@632 cba15fb0-1239-40c8-b417-11db7ca47a34

src/SampleHandler.C

@@ -3,7 +3,7 @@
 // SampleHandler.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: SampleHandler.C,v 1.8 2007-10-19 17:37:06 fgygi Exp $
+// $Id: SampleHandler.C,v 1.9 2008-06-18 03:40:53 fgygi Exp $
 
 #if USE_XERCES
 
@@ -19,9 +19,11 @@ using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
 SampleHandler::SampleHandler(Sample& s, DoubleMatrix& gfdata,
+  int& nx, int& ny, int& nz,
   vector<vector<vector<double> > > &dmat,
   Wavefunction& wfvtmp) :
-  s_(s), gfdata_(gfdata), dmat_(dmat), read_wf(false), read_wfv(false),
+  s_(s), gfdata_(gfdata), nx_(nx), ny_(ny), nz_(nz),
+  dmat_(dmat), read_wf(false), read_wfv(false),
   wfvtmp_(wfvtmp) {}
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -61,12 +63,12 @@ StructureHandler* SampleHandler::startSubHandler(const XMLCh* const uri,
   else if ( qnm == "wavefunction" )
   {
     read_wf = true;
-    return new WavefunctionHandler(s_.wf,gfdata_,dmat_);
+    return new WavefunctionHandler(s_.wf,gfdata_,nx_,ny_,nz_,dmat_);
   }
   else if ( qnm == "wavefunction_velocity" )
   {
     read_wfv = true;
-    return new WavefunctionHandler(wfvtmp_,gfdata_,dmat_);
+    return new WavefunctionHandler(wfvtmp_,gfdata_,nx_,ny_,nz_,dmat_);
   }
   else
   {

src/SampleHandler.h

@@ -3,7 +3,7 @@
 // SampleHandler.h
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: SampleHandler.h,v 1.7 2007-10-19 17:37:06 fgygi Exp $
+// $Id: SampleHandler.h,v 1.8 2008-06-18 03:40:53 fgygi Exp $
 
 #ifndef SAMPLEHANDLER_H
 #define SAMPLEHANDLER_H
@@ -27,6 +27,9 @@ class SampleHandler : public StructureHandler
   public:
 
   bool read_wf,read_wfv;
+  int& nx_;
+  int& ny_;
+  int& nz_;
 
   // Start of the root element in the structure being handled
   virtual void startElement(const XMLCh* const uri,const XMLCh* const localname,
@@ -46,7 +49,7 @@ class SampleHandler : public StructureHandler
     const XMLCh* const localname, const XMLCh* const qname,
     const StructureHandler* const subHandler);
 
-  SampleHandler(Sample& s, DoubleMatrix& gfdata,
+  SampleHandler(Sample& s, DoubleMatrix& gfdata, int& nx, int& ny, int& nz,
                 std::vector<std::vector<std::vector<double> > > &dmat,
                 Wavefunction& wfvtmp);
   ~SampleHandler();

src/SampleReader.C

@@ -3,7 +3,7 @@
 // SampleReader.C:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: SampleReader.C,v 1.27 2008-06-06 00:10:58 fgygi Exp $
+// $Id: SampleReader.C,v 1.28 2008-06-18 03:40:53 fgygi Exp $
 
 
 #include "Sample.h"
@@ -61,6 +61,7 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
   Wavefunction* wfvtmp = new Wavefunction(s.wf);
   Wavefunction* current_wf = &s.wf;
   vector<vector<vector<double> > > dmat;
+  int nx, ny, nz; // size of <grid> in wavefunction
   const int nspin = 1;
   const int current_ispin = 0;
   int current_ikp;
@@ -161,7 +162,8 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
     parser->setFeature(XMLUni::fgSAX2CoreNameSpacePrefixes, namespacePrefixes);
 
     int errorCount = 0;
-    SampleHandler* s_handler = new SampleHandler(s,gfdata,dmat,*wfvtmp);
+    SampleHandler* s_handler =
+      new SampleHandler(s,gfdata,nx,ny,nz,dmat,*wfvtmp);
 
     try
     {
@@ -208,6 +210,9 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
     if ( read_wfv )
       cout << " wavefunction velocity was read" << endl;
 
+    cout << " SampleReader::readSample: grid nx,ny,nz="
+         << nx << " " << ny << " " << nz << endl;
+
     delete s_handler;
     delete parser;
     XMLPlatformUtils::Terminate();
@@ -294,6 +299,12 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
         UnitCell uc(a,b,c);
 
         // grid
+        int ibuf[3];
+        ctxt_.ibcast_recv(3,1,ibuf,3,0,0);
+        nx = ibuf[0];
+        ny = ibuf[1];
+        nz = ibuf[2];
+
         // receive only computed ecut
         double ecut;
         ctxt_.dbcast_recv(1,1,&ecut,1,0,0);
@@ -337,6 +348,12 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
         UnitCell uc(a,b,c);
 
         // grid
+        int ibuf[3];
+        ctxt_.ibcast_recv(3,1,ibuf,3,0,0);
+        assert(nx == ibuf[0]);
+        assert(ny == ibuf[1]);
+        assert(nz == ibuf[2]);
+
         // receive only computed ecut
         double ecut;
         ctxt_.dbcast_recv(1,1,&ecut,1,0,0);
@@ -371,7 +388,7 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
           // receive grid_functions
           SlaterDet* sd = current_wf->sd(current_ispin,current_ikp);
           const Basis& basis = sd->basis();
-          FourierTransform ft(basis,basis.np(0),basis.np(1),basis.np(2));
+          FourierTransform ft(basis,nx,ny,nz);
           const int wftmpr_size = basis.real() ? ft.np012loc() :
             2*ft.np012loc();
           valarray<double> wftmpr(wftmpr_size);
@@ -456,7 +473,7 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
           // copy density matrix information
           sd->set_occ(dmat[ispin][ikp]);
           const Basis& basis = sd->basis();
-          FourierTransform ft(basis,basis.np(0),basis.np(1),basis.np(2));
+          FourierTransform ft(basis,nx,ny,nz);
 #if DEBUG
           cout << ctxt_.mype() << ": ft.np012loc()=" << ft.np012loc() << endl;
           cout << ctxt_.mype() << ": ft.context(): " << ft.context();
@@ -535,7 +552,7 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
         {
           SlaterDet* sd = wfvtmp->sd(ispin,ikp);
           const Basis& basis = sd->basis();
-          FourierTransform ft(basis,basis.np(0),basis.np(1),basis.np(2));
+          FourierTransform ft(basis,nx,ny,nz);
           //cout << ctxt_.mype() << ": ft.np012loc()=" << ft.np012loc() << endl;
           //cout << ctxt_.mype() << ": ft.context(): " << ft.context();
 

src/WavefunctionHandler.C

@@ -3,7 +3,7 @@
 // WavefunctionHandler.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: WavefunctionHandler.C,v 1.16 2008-02-12 05:39:18 fgygi Exp $
+// $Id: WavefunctionHandler.C,v 1.17 2008-06-18 03:40:53 fgygi Exp $
 
 #if USE_XERCES
 
@@ -25,9 +25,9 @@ using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
 WavefunctionHandler::WavefunctionHandler(Wavefunction& wf,
-  DoubleMatrix& gfdata,
+  DoubleMatrix& gfdata, int& nx, int& ny, int& nz,
   vector<vector<vector<double> > > &dmat) :
-  wf_(wf), gfdata_(gfdata), dmat_(dmat), ecut(0.0)
+  wf_(wf), gfdata_(gfdata), dmat_(dmat), ecut(0.0), nx_(nx), ny_(ny), nz_(nz)
 {
   // if the data is read from a file, gfdata has a finite size
   // since the grid functions were processed by XMLGFPreprocessor
@@ -177,7 +177,6 @@ void WavefunctionHandler::startElement(const XMLCh* const uri,
     //cout << " WavefunctionHandler::startElement: reference_domain" << endl;
     ruc.set(a,b,c);
     //cout << ruc;
-
   }
   else if ( locname == "density_matrix")
   {
@@ -215,18 +214,25 @@ void WavefunctionHandler::startElement(const XMLCh* const uri,
       istringstream stst(attrval);
       if ( attrname == "nx")
       {
-        stst >> nx;
+        stst >> nx_;
       }
       else if ( attrname == "ny" )
       {
-        stst >> ny;
+        stst >> ny_;
       }
       else if ( attrname == "nz" )
       {
-        stst >> nz;
+        stst >> nz_;
       }
     }
 
+    // notify listening nodes
+    int ibuf[3];
+    ibuf[0] = nx_;
+    ibuf[1] = ny_;
+    ibuf[2] = nz_;
+    wf_.context().ibcast_send(3,1,ibuf,3);
+
     if ( ecut == 0.0 )
     {
       // ecut attribute was not specified. Infer from grid size
@@ -235,9 +241,9 @@ void WavefunctionHandler::startElement(const XMLCh* const uri,
       // 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 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;
@@ -245,7 +251,6 @@ void WavefunctionHandler::startElement(const XMLCh* const uri,
       ecut = max(max(ecut0,ecut1),ecut2);
       cout << " ecut=" << 2*ecut << " Ry" << endl;
     }
-
     // notify listening nodes of ecut
     wf_.context().dbcast_send(1,1,&ecut,1);
 
@@ -298,13 +303,7 @@ void WavefunctionHandler::startElement(const XMLCh* const uri,
     assert(current_size==wf_.sd(current_ispin,current_ikp)->nst());
 
     const Basis& basis = wf_.sd(current_ispin,current_ikp)->basis();
-
-    // 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));
+    ft = new FourierTransform(basis,nx_,ny_,nz_);
     wftmp.resize((ft->np012loc()));
   }
   else if ( locname == "grid_function")
@@ -360,10 +359,11 @@ void WavefunctionHandler::endElement(const XMLCh* const uri,
   }
   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());
+    // 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_ );
 
     if ( read_from_gfdata )
     {

src/WavefunctionHandler.h

@@ -3,7 +3,7 @@
 // WavefunctionHandler.h
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: WavefunctionHandler.h,v 1.10 2007-10-19 17:37:06 fgygi Exp $
+// $Id: WavefunctionHandler.h,v 1.11 2008-06-18 03:40:53 fgygi Exp $
 
 #ifndef WavefunctionHANDLER_H
 #define WavefunctionHANDLER_H
@@ -26,7 +26,9 @@ class WavefunctionHandler : public StructureHandler
   double ecut;
   // dmat[ispin][ikp][i]
   std::vector<std::vector<std::vector<double> > > &dmat_;
-  int nx,ny,nz;
+  int& nx_;
+  int& ny_;
+  int& nz_;
   int current_gf_nx,current_gf_ny,current_gf_nz;
   std::string current_gf_encoding;
   int current_ispin,current_ikp,current_n,current_igf;
@@ -60,6 +62,7 @@ class WavefunctionHandler : public StructureHandler
     const StructureHandler* const subHandler);
 
   WavefunctionHandler(Wavefunction& wf, DoubleMatrix& gfdata,
+    int& nx, int& ny, int& nz,
     std::vector<std::vector<std::vector<double> > > &dmat);
   ~WavefunctionHandler();
 };