QCModelMatrices.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 <sstream>
00020 
00021 #include "QCCommon.hpp"
00022 #include "QCModelMatrices.hpp"
00023 #include "QCGlobalSystem.hpp"
00024 #include "QCMDSystem.hpp"
00025 #include "QCDistMDSystem.hpp"
00026 #include "QCMndo.hpp"
00027 #include "QCMndoParam.hpp"
00028 #include "QCAm1.hpp"
00029 #include "QCAm1Param.hpp"
00030 #include "QCPm3.hpp"
00031 #include "QCPm3Param.hpp"
00032 #include "QCSCFAlgo.hpp"
00033 #include "QCDCAlgo.hpp"
00034 #include "QCManager.hpp"
00035 #include "QCIntgWriter.hpp"
00036 #include "QCIntgReader.hpp"
00037 
00038 static const char* IdempotFromEnv = getenv("IDEMPOT");
00039 
00040 QCModelMatrices::QCModelMatrices (void) :
00041   overlapS(NULL),
00042   hamiltonH(),
00043   densityP(),
00044   densityPnm1(),
00045   fockF(),
00046   fockFTild(NULL),
00047   densityPTild(NULL),
00048   eigenVectCt(),
00049   eigenVal(NULL),
00050   degener(NULL),
00051   weightedDensityP(NULL),
00052   squareDensityP(NULL),
00053   threeIdMDensityP(NULL),
00054   repInterElecIntegrals(NULL),
00055   dRepInterElecIntgOnDCoord(NULL),
00056   dOverlapSOnDr(NULL) 
00057 {
00058 
00064   overlapS  = new QCOverlapGtr<QCSymMatrix>();
00065 
00066 }
00067 
00071 QCModelMatrices::~QCModelMatrices (void) {
00072 
00073   // on n'oublie jamais de detruire les alloc dyn.
00074   if (overlapS) { 
00075     delete overlapS; overlapS = NULL; 
00076   }
00077   if (eigenVal) { 
00078     delete [] eigenVal;   eigenVal = NULL; 
00079   }
00080   if (degener) { 
00081     delete [] degener;   degener = NULL; 
00082   }
00083   if (weightedDensityP) { 
00084     delete weightedDensityP;    weightedDensityP = NULL; 
00085   }  
00086   if (squareDensityP) { 
00087     delete squareDensityP;   squareDensityP = NULL; 
00088   }
00089   if (threeIdMDensityP) { 
00090     delete threeIdMDensityP; threeIdMDensityP = NULL; 
00091   }
00092 
00093   /*************************/
00094   if (fockFTild) { 
00095     delete fockFTild;   fockFTild = NULL; 
00096   }  
00097   if (densityPTild) { 
00098     delete densityPTild;  densityPTild = NULL; 
00099   }
00100   /*************************/
00101   
00102   if (repInterElecIntegrals) { 
00103     delete repInterElecIntegrals;  repInterElecIntegrals = NULL; 
00104   }
00105   if (dRepInterElecIntgOnDCoord) { 
00106     delete dRepInterElecIntgOnDCoord;   dRepInterElecIntgOnDCoord = NULL; 
00107   }
00108   if (dOverlapSOnDr) { 
00109     delete dOverlapSOnDr; dOverlapSOnDr = NULL; 
00110   }
00111 }
00112 
00116 template <class TPManager>
00117 void
00118 QCModelMatrices::allocateMatrices (TPManager& manager, int nbAOofSystem) {
00119 
00120   typedef typename TPManager::TParam TParam;
00121   
00122   const QCGeneralData& data = manager.getGeneralData();
00123 
00124   bool useOptDump = ( data.getMainAlgorithm() == QC_OPTIMAL_DAMPING ||
00125                       data.getMainAlgorithm() == QC_DC_OPTIMAL_DAMPING );
00126   bool isDCComputation = manager.isDCComputation();
00127   bool useDerivation = data.getComputationType() == QC_ENERGY_DERIV;
00128   int  nbAOspOfAtom = ( TParam::NB_S_ORB + TParam::NB_P_ORB );
00129                        
00130 
00131   overlapS->sizeMatrix(nbAOofSystem);
00132   hamiltonH.allocateMatrix(nbAOofSystem);
00133   densityP.allocateMatrix(nbAOofSystem);
00134 
00135 #ifdef DEV_SCF_CLASS_NO_MEMORY
00136   densityPnm1.allocateMatrix(nbAOofSystem);
00137 #endif
00138   fockF.allocateMatrix(nbAOofSystem);
00139   eigenVectCt.allocateMatrix(nbAOofSystem);
00140 
00141   // Si on utilise l'Optimal dumping,
00142   if (useOptDump) {
00143 
00144 #ifdef DEV_SCF_CLASS_NO_MEMORY
00145     // Il faut une matrice F_tilde et P_tilde.
00146     fockFTild    = new QCFockGtr<QCSymMatrix>(nbAOofSystem);
00147     densityPTild = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00148 #else
00149     // Les deux Matrices seront reelement allouees par
00150     // le gestionnaire de memoire
00151     fockFTild    = new QCFockGtr<QCSymMatrix>();
00152     densityPTild = new QCDensityGtr<QCSymMatrix>();
00153 #endif 
00154   
00155   }
00156   eigenVal = new QCFloat [nbAOofSystem];
00157   degener  = new int     [nbAOofSystem];
00158 
00159   // Si on fait du Divide And Conquer.
00160   if (isDCComputation) {
00161     //orbitalOccupN = new QCFloat[nbAOofSystem];
00162 
00163     // La matrice densite ponderee.
00164     weightedDensityP = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00165 
00166     if (IdempotFromEnv) {
00167       // Il faut un peu de place pour le calcul de purification.
00168       squareDensityP   = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00169       threeIdMDensityP = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00170     }
00171   }
00172 
00173   repInterElecIntegrals = new QCRepInterElec (nbAOspOfAtom*nbAOspOfAtom);
00174     
00175   if (useDerivation) {
00176     // et de la matrice pour faire le calcul des 
00177     // derivees numeriques.
00178     if (data.getDerivationType() == QC_NUMERICAL_DERIVATION) {
00179       dRepInterElecIntgOnDCoord = new QCRepInterElec(nbAOspOfAtom * nbAOspOfAtom);
00180       dOverlapSOnDr             = new QCOverlapGtr<QCSymMatrix>();
00181       dOverlapSOnDr->sizeMatrix(nbAOofSystem);
00182     }
00183   }
00184 }
00185 
00189 template <class TPManager>
00190 void 
00191 QCModelMatrices::buildMatrices (TPManager& manager, typename TPManager::TSystem::QCIterator& workingSystem) {
00192 
00193   QCIntgWriter *intgWriter = repInterElecIntegrals->getIntgWriter();
00194 
00195   // Les parametres de l'application
00196   const QCGeneralData& data = manager.getGeneralData();
00197   bool isDCComputation      = manager.isDCComputation();
00198   //
00199   // Le fichier temporaire pour les integrales de repulsion.
00200   // Ce fichier est propre a chaque SD dans le cas de D&C
00201   ostringstream idstrm;
00202   if (isDCComputation) { 
00203     idstrm << workingSystem->getId(); 
00204   }
00205 
00206   string integralsTmpFile = ( QCCommon::IntegralsTmpFile + "_" +
00207                               idstrm.str() + "_" + OUTPUT_SUFFIX );
00208   //
00209   // On ouvre le fichier de stockage des integrales
00210   //    si on est en methode d acquisition indirecte
00211   //
00212   if (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE) {
00217     if (!intgWriter) {
00218       intgWriter = repInterElecIntegrals->allocIntgWriter();
00219       repInterElecIntegrals->setIntgTmpFile(integralsTmpFile);
00220     }
00221     if (intgWriter) {
00222       intgWriter->openFile (integralsTmpFile.c_str(), "wb");
00223     
00224     } else {
00225       cerr << __FILE__ << ": buildMatrices: error openning integrals file"   << endl;
00226       exit(EXIT_FAILURE);
00227     }
00228   }
00229   // Le gestionnaire de memoire
00230   QCMemory& QCRestrict memory = manager.getMemory();
00231   //
00232   // La matrice de densite
00233   //
00234   if (data.getInitDensityMatrix() == QC_DENSITY_DIAG) {
00235     densityP.initDiag(manager, workingSystem);}
00236   //  
00237   // On genere la matrice de recouvrement.
00238   //
00239   overlapS->getMatrix().takeTmpMem(memory);
00240   overlapS->computeElems(manager, workingSystem);
00241   // Evaluation de l'utilisation memoire.
00242   memory.evaluateUsedMemory ();
00243   //
00244   // On genere H, attention, il faut cleaner d'abord car le computeElems
00245   //                ne balaie pas l'integralite de la matrice
00246   hamiltonH.init();
00247   hamiltonH.computeElems(manager, workingSystem);
00248   //
00249   overlapS->getMatrix().giveBackTmpMem(memory);
00250   //
00251   if (!manager.isDCComputation()) { // ????????????
00252     // Copie de H dans F
00253     fockF.copyMatrix(hamiltonH);
00254     fockF.completeElems(manager, workingSystem, true);
00255   }
00256 }
00257 
00261 template <class TPManager>
00262 void QCModelMatrices::prepareIterations (TPManager& manager) {
00263 
00264   QCIntgReader *intgReader = repInterElecIntegrals->getIntgReader();
00265   QCIntgWriter *intgWriter = repInterElecIntegrals->getIntgWriter();
00266 
00267   // Les parametres de l'application
00268   const QCGeneralData& data = manager.getGeneralData();
00269   
00270   // Le gestionnaire de memoire
00271   QCMemory& memory = manager.getMemory();
00272 
00273   // On update eventuellement la taille du bloc de memoire temporaire
00274   memory.updateBlockTmpMem();
00275   
00276   // On ouvre le fichier de stockage des integrales
00277   // si on est en methode d acquisition indirecte.
00278   if (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE) {
00279     //
00280     // Attention, si on ne ferme pas le fichier en ecriture, alors impossible
00281     // d'acceder aux valeurs en lecture.
00282     intgWriter->closeFile();
00283     //
00284     // Attention, oblige d attendre que le fichier soit ecrit avant de
00285     // l'ouvrir en ecriture, sinon pas de mise a jour de ce qu on a dans le
00286     // fichier et on lit un fichier vide. 
00287     // Rem : On doit pouvoir pratiquer la mise a jour
00288     // d'un fichier deja ouvert via freopen si necessaire
00289     if (!intgReader) {
00290       intgReader = repInterElecIntegrals->allocIntgReader();
00291     }
00292     
00293     if (intgReader) {
00294       intgReader->openFile(repInterElecIntegrals->getIntgTmpFile().c_str(), "rb");
00295     
00296     } else {
00297       cerr << __FILE__ << ": prepareIterations: error openning integrals file" 
00298            << endl;
00299       exit(EXIT_FAILURE);
00300     }
00301   }
00302 }
00303 
00307 QCMANAGER_METH_EXPL_INST_PARAM(void QCModelMatrices::allocateMatrices,
00308                                int);
00309 
00310 QCMANAGER_ITER_METH_EXPL_INST(void QCModelMatrices::buildMatrices);
00311 
00312 QCMANAGER_METH_EXPL_INST(void QCModelMatrices::prepareIterations);

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