QCGlobalSystem.cpp

Go to the documentation of this file.
00001 //*****************************************************************************//
00002 //                                                                             //
00003 //   Copyright (c) 2001                                                        //
00004 //      INRIA                                                                  //
00005 //      54600 VILLERS LES NANCY                                                //
00006 //      France                                                                 //
00007 //                                                                             //
00008 //*****************************************************************************//
00009 //                                                                             //
00010 //               *** NOTICE OF PROPRIETARY INFORMATION ***                     //
00011 //                                                                             //
00012 // The information contained in this file is considered proprietary and the    //
00013 // exclusive property of  INRIA. This information may not be disclosed,        //
00014 // duplicated or used, in whole or in part, for  any purpose  whatsoever       //
00015 // without express written authorization from INRIA                            //
00016 //                                                                             //
00017 //*****************************************************************************//
00018 
00019 #include "QCCommon.hpp"
00020 #include "QCReader.hpp"
00021 #include "QCMndo.hpp"
00022 #include "QCAm1.hpp"
00023 #include "QCPm3.hpp"
00024 #include "QCMndoParam.hpp"
00025 #include "QCAm1Param.hpp"
00026 #include "QCPm3Param.hpp"
00027 #include "QCSCFAlgo.hpp" 
00028 #include "QCManager.hpp"
00029 #include "QCGlobalSystem.hpp"
00030 #include "QCMemory.hpp"
00031 #include "QCSymMatrix.hpp"
00032 
00033 #include "QCTools.hpp"
00034 
00035 
00036 
00037 
00038 #ifndef QC_NO_DEBUG
00039 #  define VERBOSE_DERIV
00040 #  define DEBUG_DERIV
00041 #endif
00042 
00043 template <class TPSolver>
00044 bool QCGlobalSystem<TPSolver>::checkData(const QCGeneralData& data){
00045   QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00046   //
00047   QC_TRACE_INIT("END   QCGlobalSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00048   return true ;
00049 }
00053 template <class TPSolver>
00054 template <class TPManager>
00055 void 
00056 QCGlobalSystem<TPSolver>::computeEnergyDerivatives (TPManager& manager) {
00057 
00058 #ifdef VERBOSE_DERIV
00059   cout << endl << " ** COMPUTE ENERGY GRADIENT ** " << endl;
00060 #endif
00061   //
00062   QCMemory& QCRestrict memory = manager.getMemory();
00063   //
00064   const QCGeneralData& data = manager.getGeneralData();
00065   //
00066   typedef typename TPManager::TModel         TModel;
00067   TModel&  QCRestrict model = manager.getModel();
00068   //
00069   QCFloat deltaQ = data.getDeltaQ();
00070   //
00071   QCFloat * QCRestrict dRepIntegralSSVector = NULL;
00072 
00073   // L'iterateur
00074   typename TPManager::TSystem::QCIterator systemIter = begin();
00075 
00076   QCPoint3D shift;
00077 
00078   dRepIntegralSSVector = memory.template takeTmpMem<QCFloat>(nbAtoms*(nbAtoms-1)/2);
00079   memset(dRepIntegralSSVector, 0, (nbAtoms*(nbAtoms-1))/2*sizeof(QCFloat));
00080   overlapS->getMatrix().takeTmpMem(memory);
00081   dOverlapSOnDr->getMatrix().takeTmpMem(memory);
00082 
00083   //
00084   // ELEC ENERGY DERIV
00085   //
00086   
00087 #ifdef DEBUG_DERIV
00088   displayDisps(cout);
00089 #endif
00090 
00091   // Derivation par rapport a x.
00092   shift[COORDX] = deltaQ;
00093   shift[COORDY] = QC_ZERO;
00094   shift[COORDZ] = QC_ZERO;
00095   dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00096   model.deriveElecEnergy(manager,systemIter, shift, dRepIntegralSSVector);
00097 
00098 
00099   // Derivation par rapport a y.
00100   shift[COORDX] = QC_ZERO;
00101   shift[COORDY] = deltaQ;
00102   shift[COORDZ] = QC_ZERO;
00103   dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00104   model.deriveElecEnergy(manager,systemIter,shift,dRepIntegralSSVector);
00105 
00106 
00107   // Derivation par rapport a z.
00108   shift[COORDX] = QC_ZERO;
00109   shift[COORDY] = QC_ZERO;
00110   shift[COORDZ] = deltaQ;
00111   dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00112   model.deriveElecEnergy(manager, systemIter,shift, dRepIntegralSSVector);
00113   //
00114   dOverlapSOnDr->getMatrix().giveBackTmpMem(memory);
00115   overlapS->getMatrix().giveBackTmpMem(memory);
00116 
00117   // On divise par 2 deltaQ uniquement la derivation de l'E electronique.
00118   divDisps(QC_TWO*deltaQ);
00119 
00120   model.deriveCoreEnergy(manager,systemIter,dRepIntegralSSVector);
00121 
00122   memory.template giveBackTmpMem<QCFloat>(dRepIntegralSSVector);
00123 
00124 #ifdef DEBUG_DERIV
00125   displayDisps(cout);
00126 #endif
00127 }
00128 
00129 template <class TPSolver>
00130 void 
00131 QCGlobalSystem<TPSolver>::writeDensityOnFileAscii(const QCFiles & files,  const string& path){
00132   QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeDensityOnFileAscii   ");
00133   //
00134   std::ofstream out ;
00135   std::string fileName;
00136   //
00137   fileName =  path + "/" + files.getResultFile() + "-density.ascii";
00138   //
00139   out.open(fileName.c_str());
00140   //
00141   out.setf(std::ios_base::scientific) ;   out.precision(20) ;
00142   //
00143   QCSymMatrix& QCRestrict  D = QCModelMatrices::densityP.getMatrix();
00144   out << "BASIC "<<  D.getDim() << "  "<< D.getNbElems() << std::endl;
00145   for (int i = 0; i < D.getDim(); ++i) {
00146     for (int j = 0; j <=i; ++j) {
00147       out << "  "<< i << "  "<< j << " "<< D[i][j] << std::endl;
00148     }
00149   }  
00150   out.close();
00151   QC_TRACE_OUT("END   QCGlobalSystem<TPSolver>::writeDensityOnFileAscii   " ) ;
00152 }
00153 template <class TPSolver>
00154 void 
00155 QCGlobalSystem<TPSolver>::writeDensityOnFileBinary(const QCFiles & files,  const string& path){
00156   QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeDensityOnFileBinary  ");
00157   std::ofstream out ;
00158   std::string fileName,type("SEQ ");
00159   fileName =  path + "/" + files.getResultFile() + "-density.bin";
00160   out.open(fileName.c_str(), ios::out | ios::binary);
00161   int val[2];
00162   //
00163 
00164   QCSymMatrix& QCRestrict  D = QCModelMatrices::densityP.getMatrix();
00165   val[0] = D.getDim(); val[1] = D.getNbElems()*sizeof(QCFloat) ;
00166   out.write (type.c_str(), 4);
00167   out.write ((char*)(&val[0]), 2*sizeof(int));
00168   out.write ((char*)(D.getElems()), val[1]);
00169 
00170   out.close();
00171   QC_TRACE_OUT("END   QCGlobalSystem<TPSolver>::writeDensityOnFileBinary " ) ;
00172 }
00173 template <class TPSolver>
00174 void 
00175 QCGlobalSystem<TPSolver>::readDensityFromFile(const QCFiles & files,  const string& path){
00176   //
00177   //
00178   QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFile   ");
00179   std::string fileName;
00180   fileName = files.getDensityFile();
00181   //
00182   //  QCSymMatrix& QCRestrict  D = QCModelMatrices::densityP.getMatrix();
00183   if(fileName.find(".bin")!= std::string::npos){
00184     this->readDensityFromFileBinary( files, path);
00185   }
00186   else if(fileName.find(".ascii",0)!= std::string::npos){
00187     readGlobalDensityFromAsciiFile(fileName,QCModelMatrices::densityP.getMatrix());
00188   }
00189   else {
00190     std::cerr <<" Wrong file for the density matrix. the extension must be .bin or .ascii"<<std::endl;
00191     exit(EXIT_FAILURE) ;
00192   }
00193   QC_TRACE_INIT("END   QCGlobalSystem<TPSolver>::readDensityFromFile   ");
00194 }
00195 
00196 template <class TPSolver>
00197 void 
00198 QCGlobalSystem<TPSolver>::readDensityFromFileAscii(const QCFiles & files,  const string& path){
00199   QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFileAscii   ") ;
00200   std::ifstream data ;
00201   std::string fileName, type;
00202   int dim = 0,  nbelement = 0;
00203   fileName = files.getDensityFile();
00204   data.open(fileName.c_str());
00205   if (!data){
00206     std::cerr << "Error to open density Matrix in File " << fileName <<std::endl; 
00207     exit(EXIT_FAILURE) ;
00208   }
00209  //
00210   QCSymMatrix& QCRestrict  D = QCModelMatrices::densityP.getMatrix();
00211   
00212   data >> type >>  dim >> nbelement;
00213   if(type != "SEQ") {
00214     std::cerr << "Bad type of file " << type << " must be seq FOR binary file "<<std::endl;
00215     exit (EXIT_FAILURE) ;
00216   }
00217   if ( dim != D.getDim() ){
00218     std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
00219               << " Dimension in File : " << dim <<"  Dimension given by the model" << D.getDim() <<std::endl ;
00220     exit (EXIT_FAILURE) ;
00221   }
00222   int i, j ;
00223   QCFloat val ;
00224   for (int k = 0; k < nbelement ; ++k) {
00225     data >> i >> j >>val ;
00226     if (i < j ){ 
00227       std::cerr <<" Error we must read a Symetric Matrix " <<std::endl ;
00228       exit(EXIT_FAILURE) ;
00229     }
00230     D[i][j] = val ;
00231   }
00232   data.close();
00233   QC_TRACE_INIT("END   QCGlobalSystem<TPSolver>::readDensityFromFileAscii   ")
00234 }
00235 template <class TPSolver>
00236 void 
00237 QCGlobalSystem<TPSolver>::readDensityFromFileBinary(const QCFiles & files,  const string& path){
00238   QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFileBinary  ") ;
00239   std::ifstream data ;
00240   std::string fileName, tt ;
00241   char type[4];
00242   fileName = files.getDensityFile();
00243   data.open(fileName.c_str(), ios::in|ios::binary);
00244   if (!data){
00245    std::cerr << "Error to open density Matrix in File " << fileName <<std::endl; 
00246     exit(EXIT_FAILURE) ;
00247   }
00248   //
00249   QCSymMatrix& QCRestrict  D = QCModelMatrices::densityP.getMatrix();
00250   int val[2] ;
00251   //
00252   data.read(type, 4); tt = type ;
00253   data.read((char*)(&val[0]), 2*sizeof(int));
00254   if(tt != "SEQ") {
00255     std::cerr << "Bad type of file " << type << " must be seq FOR binary file "<<std::endl;
00256     exit (EXIT_FAILURE) ;
00257   }
00258   if ( val[0] != D.getDim() ){
00259     std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
00260               << " Dimension in File : " << val[0] <<"  Dimension given by the model" << D.getDim() <<std::endl ;
00261     exit (EXIT_FAILURE) ;
00262   }
00263   data.read((char*)(D.getElems()), val[1]);
00264   //
00265   data.close();
00266   QC_TRACE_INIT("END   QCGlobalSystem<TPSolver>::readDensityFromFileBinary   ")
00267 }
00268 template <class TPSolver>
00269 void 
00270 QCGlobalSystem<TPSolver>::writeSystemOnFile(const QCFiles & files,  const string& path){
00271   QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeSystemOnFile   ") ;
00272   std::ofstream out ;
00273   std::string fileName;
00274   //
00275   fileName = path + "/" + files.getResultFile() + ".system";
00276   out.open(fileName.c_str());
00277   //
00278   QCSystem::writeQCIAtomicSystem(out) ;
00279   QC_TRACE_OUT("END   QCGlobalSystem<TPSolver>::writeSystemOnFile  ") ;
00280     }
00284 template class QCGlobalSystem<QCSCFAlgo>;
00285 
00286 
00287 QCMANAGER_METH_EXPL_INST_GLOBAL(void QCGlobalSystem<QCSCFAlgo>::computeEnergyDerivatives);

Generated on Sat Jan 28 21:07:26 2006 for QC++ by  doxygen 1.4.4