QMcBeaver: QMCConfigIO.cpp Source File

00001 /*
00002   Copyright (c) Amos G. Anderson 2005
00003   Distributed under GNU general public license (GPL)
00004   No guarantee or warantee regarding usability or stability is expressed or implied.
00005   nitroamos@gmail.com
00006 */
00007 
00008 #include "QMCConfigIO.h"
00009 
00026 static const bool inBinary = true;
00027 
00028 int QMCConfigIO::numRead      = 0;
00029 int QMCConfigIO::numWritten   = 0;
00030 int QMCConfigIO::numElectrons = 0;
00031 bool QMCConfigIO::areWriting  = false;
00032 string QMCConfigIO::filename  = "";
00033 
00034 #ifdef USEHDF5
00035 H5File * QMCConfigIO::h5_f      = 0;
00036 DataSet * QMCConfigIO::dst_r    = 0;
00037 DataSet * QMCConfigIO::dst_g    = 0;
00038 DataSet * QMCConfigIO::dst_o    = 0;
00039 CompType * QMCConfigIO::ct      = 0;
00040 #else
00041 fstream * QMCConfigIO::config_strm = 0;
00042 #endif
00043 
00044 QMCConfigIO::QMCConfigIO()
00045 {
00046 #ifdef USEHDF5
00047   h5_f  = 0;
00048   dst_r = 0;
00049   dst_o = 0;
00050   dst_g = 0;
00051   ct = 0;
00052 #else
00053   config_strm = 0;
00054 #endif
00055   areWriting = true;
00056   numWritten = 0;
00057   numRead    = 0;
00058 }
00059 
00060 QMCConfigIO::~QMCConfigIO()
00061 {
00062   close();
00063 }
00064 
00065 QMCConfigIO::QMCConfigIO(int numE)
00066 {
00067 #ifdef USEHDF5
00068   h5_f  = 0;
00069   dst_r = 0;
00070   dst_o = 0;  
00071   dst_g = 0;  
00072   ct = 0;
00073 #else
00074   if(!inBinary)
00075     {
00076       clog << "Warning: CFGS file is not binary, severe performance penalty!!" << endl;
00077     }
00078 
00079   config_strm = 0;
00080 #endif
00081   areWriting = true;
00082   numWritten = 0;
00083   numRead    = 0;
00084   filename   = string("");
00085   numElectrons = numE;
00086 }
00087 
00088 void QMCConfigIO::open(string name, bool forWriting)
00089 {
00090 #ifdef USEHDF5
00091   //I think it's a good idea to use the common HDF5 file extension
00092   name += ".h5";
00093   if( name == filename && areWriting == forWriting && h5_f != 0 )
00094     return;
00095 #else
00096   if( name == filename && areWriting == forWriting && config_strm != 0 )
00097     return;
00098 #endif
00099 
00100   areWriting = forWriting;
00101   filename = name;
00102   close();
00103 
00104   if(forWriting)
00105     numWritten = 0;
00106 
00107   if(filename == "")
00108   {
00109     clog << "Warning: empty filename" << endl;
00110     return;
00111   }
00112   
00113 #ifdef USEHDF5
00114 
00115   if( !forWriting && !H5File::isHdf5(filename) )
00116     {
00117       clog << "Error: the file " << filename << " is not an HDF5 file!" << endl;
00118     }
00119 
00120   bool printError = false;
00121   try {
00122     //Exception::dontPrint();
00123 
00124     ct = new CompType(sizeof(ct_typ));
00125     ct->insertMember("SCF_Lap",  HOFFSET(ct_typ, SCF_Lap ), PredType::NATIVE_DOUBLE);
00126     ct->insertMember("weight",   HOFFSET(ct_typ, weight  ), PredType::NATIVE_DOUBLE);
00127     ct->insertMember("PE",       HOFFSET(ct_typ, PE      ), PredType::NATIVE_DOUBLE);
00128     ct->insertMember("lnJ",      HOFFSET(ct_typ, lnJ     ), PredType::NATIVE_DOUBLE);
00129 
00130     const hsize_t dim = 3*numElectrons;
00131 
00132     /*
00133       See User Guide, Chapter 2, section 8 for descriptions
00134       the "core" driver keeps it all in RAM, paramter is how much memory to increment
00135       the "sec2" driver is default, minimal buffering      
00136     */
00137     FileAccPropList   access;
00138     //access.setCore(1024, (hbool_t)false);
00139     //access.setStdio();
00140 
00141     FileCreatPropList create;
00142 
00143     if(forWriting)
00144       {
00145         h5_f = new H5File(filename, H5F_ACC_TRUNC,  create.getId(), access.getId());
00146         DSetCreatPropList dcpl_r, dcpl_o, dcpl_g;
00147 
00148         /*
00149           Higher number corresponds to more compression.
00150           0 through 9 allowed
00151         */
00152         const int compressionLevel = 0;
00153         
00154         /*
00155           The larger this parameter is, the more HDF5 will buffer
00156           before writting. Testing on matrix suggests this only
00157           makes a couple % of difference.
00158         */
00159         const hsize_t chunkFactor = 100;
00160 
00161 
00162         //what's the best way to order the dimensions?
00163         //which one should be unlimited?
00164         hsize_t max_rgr_size[2] = {H5S_UNLIMITED, dim};
00165         hsize_t     rgr_size[2] = {chunkFactor,   dim};
00166         DataSpace dsp1(2,rgr_size,max_rgr_size);
00167         dcpl_r.setChunk(2,rgr_size);
00168         dcpl_r.setDeflate(compressionLevel);
00169         dcpl_r.setFillTime(H5D_FILL_TIME_NEVER);
00170         dcpl_r.setAllocTime(H5D_ALLOC_TIME_LATE);
00171         dst_r = new DataSet(h5_f->createDataSet("R",PredType::NATIVE_DOUBLE,dsp1,dcpl_r));
00172         dsp1.close();
00173 
00174         DataSpace dsp3(2,rgr_size,max_rgr_size);
00175         dcpl_g.setChunk(2,rgr_size);
00176         dcpl_g.setDeflate(compressionLevel);
00177         dcpl_g.setFillTime(H5D_FILL_TIME_NEVER);
00178         dcpl_g.setAllocTime(H5D_ALLOC_TIME_LATE);
00179         dst_g = new DataSet(h5_f->createDataSet("Grad",PredType::NATIVE_DOUBLE,dsp3,dcpl_g));
00180         dsp3.close();
00181         
00182         hsize_t max_oth_size[1] = {H5S_UNLIMITED};
00183         hsize_t     oth_size[1] = {chunkFactor};
00184         DataSpace dsp2(1,oth_size,max_oth_size);
00185         dcpl_o.setChunk(1,oth_size);
00186         dcpl_o.setDeflate(compressionLevel);
00187         dcpl_o.setFillTime(H5D_FILL_TIME_NEVER);
00188         dcpl_o.setAllocTime(H5D_ALLOC_TIME_LATE);
00189         dst_o = new DataSet(h5_f->createDataSet("Others",*ct,dsp2,dcpl_o)); 
00190         dsp2.close();
00191 
00192         numWritten = 0;
00193 
00194       } else {
00195         h5_f = new H5File(filename, H5F_ACC_RDONLY, create.getId(), access.getId());
00196         dst_r = new DataSet(h5_f->openDataSet("R")); 
00197         dst_g = new DataSet(h5_f->openDataSet("Grad")); 
00198         dst_o = new DataSet(h5_f->openDataSet("Others")); 
00199         numRead = 0;
00200       }
00201   }
00202 
00203   // catch failure caused by the H5File operations
00204   catch( FileIException error )
00205     {
00206       error.printError();
00207       printError = true;
00208     }
00209   
00210   // catch failure caused by the DataSet operations
00211   catch( DataSetIException error )
00212     {
00213       error.printError();
00214       printError = true;
00215     }
00216   
00217   // catch failure caused by the DataSpace operations
00218   catch( DataSpaceIException error )
00219     {
00220       error.printError();
00221       printError = true;
00222     }
00223   
00224   // catch failure caused by the DataSpace operations
00225   catch( DataTypeIException error )
00226     {
00227       error.printError();
00228       printError = true;
00229     }
00230 
00231   if(printError)
00232     {
00233       cout << "Error: HDF5 can't open file " << filename << endl;
00234     }
00235 
00236 #else  
00237   ios_base::openmode mode;
00238   
00239   if(forWriting)
00240     {
00241       mode = ios_base::trunc;
00242       mode |= ios_base::out;
00243     }
00244   else 
00245     {
00246       mode = ios_base::in;
00247     }
00248   
00249   if(inBinary) mode |= ios_base::binary;
00250   
00251   config_strm = new fstream(filename.c_str(), mode);
00252   
00253   if( !(*config_strm) || config_strm->bad() )
00254     {
00255       cerr << "Error: Can't open " << filename;
00256       if(config_strm->rdstate() & ios::badbit)
00257         cerr << " badbit ";
00258       if(config_strm->rdstate() & ios::failbit)
00259         cerr << " failbit ";
00260       if(config_strm->rdstate() & ios::eofbit)
00261         cerr << " eofbit ";
00262       cerr << endl;
00263       exit(1);
00264     }
00265   
00266   //this is just to make sure that the EOF flag is cleared
00267   config_strm->clear();
00268 #endif
00269 }
00270 
00271 void QMCConfigIO::writeCorrelatedSamplingConfiguration(Array2D<double> & R,
00272                                                        double SCF_Laplacian_PsiRatio,
00273                                                        Array2D<double> & SCF_Grad_PsiRatio,
00274                                                        double lnJastrow,
00275                                                        double PE,
00276                                                        double weight)
00277 {
00278   if(!areWriting)
00279     {
00280       cerr << "Warning: the file is open for reading, not writing." << endl;
00281       return;
00282     }
00283 
00284 #ifdef USEHDF5
00285 
00286   packet.SCF_Lap = SCF_Laplacian_PsiRatio;
00287   packet.weight  = weight;
00288   packet.PE      = PE;
00289   packet.lnJ     = lnJastrow;
00290   
00291   hsize_t      dim       = 3*numElectrons;
00292   hsize_t      dims1[2]  = {1,              dim};  
00293   hsize_t      size1[2]  = {numWritten+100, dim};
00294   hsize_t      offs1[2]  = {numWritten,       0};
00295 
00296   dst_r->extend( size1 );
00297   DataSpace ms1(2, dims1);
00298   DataSpace fs1,fs3;
00299   fs1 = dst_r->getSpace();
00300   fs1.selectHyperslab(H5S_SELECT_SET, dims1, offs1);
00301   dst_r->write(R.array(), PredType::NATIVE_DOUBLE, ms1, fs1);
00302 
00303   dst_g->extend( size1 );
00304   DataSpace ms3(2, dims1);
00305   fs3 = dst_g->getSpace();
00306   fs3.selectHyperslab(H5S_SELECT_SET, dims1, offs1);
00307   dst_g->write(SCF_Grad_PsiRatio.array(), PredType::NATIVE_DOUBLE, ms3, fs3);
00308   fs3.close();
00309   ms3.close();
00310 
00311   hsize_t      dims2[1]  = {1};  
00312   hsize_t      size2[1]  = {numWritten+100};
00313   hsize_t      offs2[1]  = {numWritten};
00314   dst_o->extend( size2 );
00315 
00316   DataSpace ms2(1, dims2);
00317   DataSpace fs2;
00318   fs2 = dst_o->getSpace();
00319   fs2.selectHyperslab(H5S_SELECT_SET, dims2, offs2);
00320   dst_o->write(&packet, (*ct), ms2, fs2);
00321   fs2.close();
00322   ms2.close();
00323   
00324 #else
00325   if(config_strm == 0)
00326   {
00327     cerr << "Warning: config_strm is zero (write sample)." << endl;
00328     return;    
00329   }
00330   
00331   if(inBinary)
00332   {
00333     config_strm->write( (char*) &numElectrons, sizeof(int) );
00334     config_strm->write( (char*) R.array(), sizeof(*R.array())*R.dim1()*3 );
00335     config_strm->write( (char*) &SCF_Laplacian_PsiRatio, sizeof(double) );
00336     config_strm->write( (char*) SCF_Grad_PsiRatio.array(),
00337                         sizeof(*SCF_Grad_PsiRatio.array())*SCF_Grad_PsiRatio.dim1()*3 );
00338     config_strm->write( (char*) &lnJastrow, sizeof(double) );
00339     config_strm->write( (char*) &PE, sizeof(double) ); 
00340     config_strm->write( (char*) &weight, sizeof(double) ); 
00341   }
00342   else
00343   {
00344     setPrecision();
00345     fstream & strm = *config_strm;
00346     
00347     strm << "&" << endl;
00348     strm << "R\t" << numElectrons << endl;
00349     for(int i=0;i<numElectrons;i++)
00350     {
00351       //now we are printing out all the electrons
00352       
00353       strm << "\t" << i << "\t";
00354       strm << R(i,0) << "\t";
00355       strm << R(i,1) << "\t";
00356       strm << R(i,2) << endl;
00357     }
00358     
00359     // Print the sum of the laplacians from the alpha and beta electrons
00360     strm << "D1" << endl;
00361     strm << "\t" << SCF_Laplacian_PsiRatio << endl;
00362     
00363     // print the GradPsiRatio excluding the Jastrow
00364     strm << "D2" << endl;
00365     
00366     for(int i=0; i<numElectrons; i++)
00367     {
00368       for(int j=0; j<3; j++)
00369       {
00370         strm << "\t" << SCF_Grad_PsiRatio(i,j);
00371       }
00372       strm << endl;
00373     }
00374     
00375     strm << "lnJ\t" << endl;
00376     strm << "\t" << lnJastrow << endl;
00377     strm << "PE\t" << endl;
00378     strm << "\t" << PE << endl;
00379     strm << "W\t" << endl;
00380     strm << "\t" << weight << endl;
00381   }
00382 #endif
00383 
00384   numWritten++;
00385 }
00386 
00387 void QMCConfigIO::readCorrelatedSamplingConfiguration(Array2D<double> & R,
00388                                                       double & SCF_Laplacian_PsiRatio,
00389                                                       Array2D<double> & SCF_Grad_PsiRatio,
00390                                                       double & lnJastrow,
00391                                                       double & PE,
00392                                                       double & weight)
00393 {
00394   if(areWriting)
00395     {
00396       cerr << "Warning: the file is open for writing, not reading." << endl;
00397       return;
00398     }
00399 
00400 #ifdef USEHDF5
00401   int err = 0;
00402 
00403   hsize_t      dim       = 3*numElectrons;
00404   hsize_t      dims1[2]  = {1,          dim};  
00405   hsize_t      offs1[2]  = {numRead,      0};
00406   offs1[0] = numRead;
00407   
00408   DataSpace ms1(2, dims1);
00409   DataSpace fs1;
00410   fs1 = dst_r->getSpace();
00411   fs1.selectHyperslab(H5S_SELECT_SET, dims1, offs1);
00412   dst_r->read(R.array(), PredType::NATIVE_DOUBLE, ms1, fs1);
00413   fs1.close();
00414 
00415   DataSpace ms3(2, dims1);
00416   DataSpace fs3;
00417   fs3 = dst_g->getSpace();
00418   fs3.selectHyperslab(H5S_SELECT_SET, dims1, offs1);
00419   dst_g->read(SCF_Grad_PsiRatio.array(), PredType::NATIVE_DOUBLE, ms3, fs3);
00420   fs3.close();
00421   
00422   hsize_t      dims2[1]  = {1};  
00423   hsize_t      offs2[1]  = {numRead};
00424   offs2[0] = numRead;
00425   
00426   DataSpace ms2(1, dims2);
00427   DataSpace fs2;
00428   fs2 = dst_o->getSpace();
00429   fs2.selectHyperslab(H5S_SELECT_SET, dims2, offs2);  
00430   dst_o->read(&packet, (*ct), ms2, fs2);
00431   fs2.close();
00432   
00433   if(err < 0) exit(0);
00434   SCF_Laplacian_PsiRatio = packet.SCF_Lap;
00435   weight                 = packet.weight;
00436   PE                     = packet.PE;
00437   lnJastrow              = packet.lnJ;
00438   
00439 #else
00440   if(config_strm == 0)
00441   {
00442     cerr << "Warning: config_strm is zero (read sample)." << endl;
00443     return;    
00444   }
00445     
00446   if(inBinary)
00447   {
00448     config_strm->read( (char*) &numElectrons, sizeof(int) );
00449     config_strm->read( (char*) R.array(), sizeof(*R.array())*R.dim1()*3 );
00450     config_strm->read( (char*) &SCF_Laplacian_PsiRatio, sizeof(double) );
00451     config_strm->read( (char*) SCF_Grad_PsiRatio.array(),
00452                         sizeof(*SCF_Grad_PsiRatio.array())*SCF_Grad_PsiRatio.dim1()*3 );
00453     config_strm->read( (char*) &lnJastrow, sizeof(double) );
00454     config_strm->read( (char*) &PE, sizeof(double) ); 
00455     config_strm->read( (char*) &weight, sizeof(double) ); 
00456   }
00457   else
00458   {
00459     string stemp;
00460     int itemp;
00461     fstream & strm = *config_strm;
00462     
00463     strm >> stemp;
00464     strm >> stemp;
00465     strm >> stemp;
00466     
00467     // Read in the electronic positions
00468     for(int i=0; i<numElectrons; i++)
00469     {
00470       strm >> itemp;
00471 
00472       strm >> R(i,0);
00473       strm >> R(i,1);
00474       strm >> R(i,2);
00475     }
00476     
00477     // read D1
00478     strm >> stemp;
00479     strm >> SCF_Laplacian_PsiRatio;
00480     
00481     // read D2
00482     strm >> stemp;
00483     for(int i=0; i<numElectrons; i++)
00484     {
00485       strm >> SCF_Grad_PsiRatio(i,0);
00486       strm >> SCF_Grad_PsiRatio(i,1);
00487       strm >> SCF_Grad_PsiRatio(i,2);
00488     }
00489     
00490     // read lnJ
00491     strm >> stemp;
00492     strm >> lnJastrow;
00493 
00494     // read PE
00495     strm >> stemp;
00496     strm >> PE;
00497 
00498     // read the weight
00499     strm >> stemp;
00500     strm >> weight;
00501   }
00502   config_strm->flush();
00503 #endif
00504 
00505   numRead++;
00506 }
00507 
00508 void QMCConfigIO::close()
00509 {
00510 #ifdef USEHDF5
00511   if(h5_f != 0){
00512     h5_f->close();
00513     delete h5_f;     h5_f  = 0;
00514     delete dst_r;     dst_r  = 0;
00515     delete dst_o;     dst_o  = 0;
00516     delete dst_g;     dst_g  = 0;
00517     delete ct;       ct    = 0;
00518   }
00519 #else
00520   if(config_strm != 0){
00521     config_strm->close();
00522     delete config_strm;
00523     config_strm = 0;
00524   }
00525 #endif
00526 }
00527 
00528 string QMCConfigIO::getFilename()
00529 {
00530   return filename;
00531 }
00532 
00533 bool QMCConfigIO::eof()
00534 {
00535 #ifdef USEHDF5
00536   int err;
00537 
00538   if(numRead < numWritten)
00539     return false;
00540   return true;
00541 
00542 #else
00543   if(config_strm == 0)
00544   {
00545     cerr << "Warning: config_strm is zero (eof)." << endl;
00546     return true;    
00547   }
00548     
00549   return config_strm->eof();
00550 #endif
00551 }
00552 
00553 void QMCConfigIO::setPrecision()
00554 {
00555 #ifdef USEHDF5
00556 #else
00557   if(config_strm == 0)
00558   {
00559     cerr << "Warning: config_strm is zero (setPrecision)." << endl;
00560     return;    
00561   }
00562   
00563   config_strm->setf(ios::fixed,ios::floatfield);
00564   config_strm->precision(10);
00565 #endif
00566 }
00567