QCFockGtr.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 
00020 #include <iomanip>
00021 
00022 #include "config.h"
00023 #include "QCCommon.hpp"
00024 #include "QCFockGtr.hpp"
00025 #include "QCSymMatrix.hpp"
00026 #include "QCMndo.hpp"
00027 #include "QCAm1.hpp"
00028 #include "QCPm3.hpp"
00029 #include "QCMndoParam.hpp"
00030 #include "QCAm1Param.hpp"
00031 #include "QCPm3Param.hpp"
00032 #include "QCGlobalSystem.hpp"
00033 #include "QCMDSystem.hpp"
00034 #include "QCDistMDSystem.hpp"
00035 #include "QCSCFAlgo.hpp"
00036 #include "QCDCAlgo.hpp"
00037 #include "QCManager.hpp"
00038 #include "QCIntgReader.hpp"
00039 
00040 
00041 
00045 template <class TPMatrix>
00046 QCFockGtr<TPMatrix>::QCFockGtr (void) : 
00047   QCMatElemGtr<TPMatrix>() 
00048 {}
00049 
00050 
00051 
00055 template <class TPMatrix>
00056 QCFockGtr<TPMatrix>::QCFockGtr (int dimMatrix) : 
00057   QCMatElemGtr<TPMatrix>(dimMatrix) 
00058 {}
00059 
00060 
00061 
00065 template <class TPMatrix>
00066 QCFockGtr<TPMatrix>::~QCFockGtr (void) {
00067 }
00068 
00069 
00070 
00071 
00079 template <class TPMatrix>
00080 template <class TPManager>
00081 void QCFockGtr<TPMatrix>::completeElems (TPManager&  manager,
00082                                          typename TPManager::TSystem::QCIterator& workingSystem,
00083                                          const bool&    isItfirstCall) {
00084 
00085   const QCGeneralData& data = manager.getGeneralData();
00086   //
00087   QCIntgAcquisitionMethod intgAcquisitionMethod = data.getIntgAcquisitionMethod();
00088   //
00089   // On recupere avec des restricts pour le calcul les
00090   // structures qui contiennent les donnees a manipuler
00091   //
00092   QCSystem&        QCRestrict currentSystem   = workingSystem->getSystem();
00093   QCModelMatrices& QCRestrict currentMatrices = workingSystem->getMatrices();
00094   // Les parametres
00095   const typename TPManager::TParam *QCRestrict parameters =    manager.getParameters();
00096   // Le modele
00097   typename TPManager::TModel& QCRestrict model = manager.getModel();
00098   // Les atomes
00099   QCPoint3D  coordsA;
00100   QCPoint3D  coordsB;
00101   // leurs parametres
00102   const typename TPManager::TParam *QCRestrict parameterA;
00103   const typename TPManager::TParam *QCRestrict parameterB;
00104   // distance entre a et b.
00105   QCFloat rAB;
00106   // La matrice densite
00107   QCSymMatrix& QCRestrict densityMatrix = currentMatrices.getDensityP().getMatrix();
00108   //
00109   // La matrice de retour de la methode de calcul des integrales de repulsion inter elec.
00110   QCRepInterElec& QCRestrict repInterElecIntegrals =  currentMatrices.getRepInterElecIntegrals();
00111   // La matrice des integrales de repulsion..
00112   QCMatrix& QCRestrict repIntegralsMatrix =  repInterElecIntegrals.getMatrix();  
00113   //
00114   // Matrice de Fock locale a cette methode. Attention,
00115   // cette matrice n'est pas symetrique.
00116   QCSymMatrix& QCRestrict fockFAA = model.getSpWorkingAA();
00117   QCMatrix&    QCRestrict fockFAB = model.getSpWorkingAB();
00118   QCSymMatrix& QCRestrict fockFBB = model.getSpWorkingBB();
00119   //
00120   // On a enfin besoin de la matrice densite
00121   QCSymMatrix& QCRestrict densityPAA = model.getSpWorkingAA2();
00122   QCSymMatrix& QCRestrict densityPBB = model.getSpWorkingBB2();
00123   QCMatrix&    QCRestrict densityPAB = model.getSpWorkingAB2();
00124   //
00125   // Les premiers index d'orbitales de A et B dans la matrice complete du
00126   // systeme et l'index du premier successeur de A tjs dans la  matrice complete.
00127   //
00128   int firstAOofA, firstAOofB;
00129   // valeur intermediaire
00130   register QCFloat interVal, interVal1, interVal2;
00131   int i, j, mu, nu, lda, sma, muNu, ldaSma;
00132   //
00133   // La double boucle sur les atomes.  je prends i < j.
00134   //
00135   for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00136     currentSystem.getPointAt(i, coordsA);
00137     parameterA = &parameters[currentSystem.getParamIndexAt(i)];
00138     firstAOofA = currentSystem.getFirstAOAt(i);
00139     //
00140     // On met le bloc de matrice qui nous interesse 
00141     // dans une sous-matrice
00142     this->matrix.extractSubTriangle(fockFAA, firstAOofA, parameterA->getNbAO());
00143     //
00144     //Idem pour la matrice densite
00145     densityMatrix.extractSubTriangle(densityPAA, firstAOofA, parameterA->getNbAO());
00146     //
00147     // On commence par le debut, soit la completions des termes diagonaux F_mu_mu.
00148     //
00149     for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00150       // Le premier terme avec une simple somme sur les nu
00151       interVal = QC_ZERO;
00152       for (nu = 0; nu < parameterA->getNbAO(); ++nu) {
00153         interVal += densityPAA[nu][nu] * (parameterA->getG(mu, nu) - QC_HALF * parameterA->getH(mu, nu) );
00154       }
00155       fockFAA[mu][mu] += interVal;
00156       for (nu = 0; nu < mu; ++nu) {
00157         fockFAA[mu][nu] += QC_HALF * densityPAA[mu][nu] * 
00158           ( QC_THREE * parameterA->getH(mu, nu) - parameterA->getG(mu, nu) );
00159       }
00160     }
00161     //
00162     // Le cas ou A est different de B. On veut assigner les F_mu_lambda
00163     // 
00164     for (j = 0; j < i; ++j) {
00165 
00166       currentSystem.getPointAt(j, coordsB);
00167       parameterB = &parameters[currentSystem.getParamIndexAt(j)];
00168       firstAOofB = currentSystem.getFirstAOAt(j);
00169       rAB        = currentSystem.distance(i,j);
00170       //
00171       // On met 2 blocs de matrices de Fock dans des sous-matrices.
00172       this->matrix.extractSubMatrix(fockFAB,firstAOofA, parameterA->getNbAO(),
00173                               firstAOofB, parameterB->getNbAO() );
00174       this->matrix.extractSubTriangle(fockFBB,firstAOofB, parameterB->getNbAO() );
00175 
00176       // On extrait aussi les 2 blocs de matrice densite.
00177       densityMatrix.extractSubMatrix(densityPAB, firstAOofA, parameterA->getNbAO(),
00178                                      firstAOofB, parameterB->getNbAO() );
00179       densityMatrix.extractSubTriangle(densityPBB, firstAOofB, parameterB->getNbAO() );
00180       //
00181       // Calcul des elements de la matrice de rep inter elec.
00182       //
00183       repInterElecIntegrals.computeElems(// Attention, le premier entree est celui dont l'indice
00184                                          // varie le plus vite
00185                                          parameterA,parameterB,coordsA,coordsB, rAB,
00186                                          // ssACalc et ssBCalc recoivent false,
00187                                          // car ici, on fait tous les calculs
00188                                          // des elements de la matrice
00189                                          false, false, false, intgAcquisitionMethod, isItfirstCall);
00190       /************
00191        * A <--> B *
00192        ************/
00193       for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00194         for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00195           interVal = QC_ZERO;
00196           muNu = mu * parameterA->getNbAO();
00197           for (nu = 0;  nu < parameterA->getNbAO(); ++nu, ++muNu) {
00198             ldaSma = lda * parameterB->getNbAO();
00199             for (sma = 0; sma < parameterB->getNbAO(); ++sma, ++ldaSma) {
00200 
00201               // Completion des termes F_mu_lda
00202               interVal += densityPAB[nu][sma] * repIntegralsMatrix[muNu][ldaSma];
00203             }
00204           }
00205           fockFAB[mu][lda] -= QC_HALF * interVal;
00206         }
00207       }
00208       /************
00209        * A <--> A *
00210        ************/
00211       
00212       
00213       if ( (data.getPartitionType() == QC_STANDARD_PART &&  workingSystem->getZone(j) != QC_SHELL2) || 
00214            (data.getPartitionType() == QC_DIXON_PART &&  workingSystem->getZone(j) == QC_CORE) ) {
00215         for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00216           muNu = mu * parameterA->getNbAO();
00217           for (nu = 0;  nu <= mu; ++nu, ++muNu) {
00218             interVal1  = QC_ZERO;
00219             interVal2  = QC_ZERO;
00220             for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00221               ldaSma = lda * parameterB->getNbAO();
00222               for (sma = 0; sma < lda; ++sma, ++ldaSma) {
00223                 
00224                 // Completion des termes F_mu_mu ou F_mu_nu de l atome A.
00225                 interVal1 += densityPBB[lda][sma] * repIntegralsMatrix[muNu][ldaSma];
00226               }
00227               // Attention, ici, je profite du fait que sma == lda en sortie
00228               // de la bcle precedente, idem pour ldaSma.
00229               interVal2 += densityPBB[lda][sma] * repIntegralsMatrix[muNu][ldaSma];
00230             }
00231             // On multiplie par 2 interVal pour avoir 2 * 
00232             // les elem sym du bloc diagonal.
00233             fockFAA[mu][nu] += QC_TWO * interVal1;
00234             fockFAA[mu][nu] += interVal2;
00235           }
00236         }
00237       }
00238       /************
00239        * B <--> B *
00240        ************/
00241 
00242       if ( (data.getPartitionType() == QC_STANDARD_PART &&  workingSystem->getZone(i) != QC_SHELL2) || 
00243            (data.getPartitionType() == QC_DIXON_PART && workingSystem->getZone(i) == QC_CORE) ) {
00244         //      
00245         for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00246           ldaSma = lda * parameterB->getNbAO();
00247           for (sma = 0; sma <= lda; ++sma, ++ldaSma) {
00248             interVal1  = QC_ZERO;
00249             interVal2  = QC_ZERO;
00250             for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00251               muNu = mu * parameterA->getNbAO();
00252               for (nu = 0;  nu < mu; ++nu, ++muNu) {
00253                 //
00254                 // Completion des termes F_mu_mu ou F_mu_nu de l atome B.             
00255                 interVal1 +=  densityPAA[mu][nu] * repIntegralsMatrix[muNu][ldaSma];
00256               }
00257               //
00258               // Attention, ici, je profite du fait que mu == nu en sortie
00259               // de la bcle precedente, idem pour muNu.
00260               interVal2 += densityPAA[mu][nu] * repIntegralsMatrix[muNu][ldaSma];
00261             }
00262             // On multiplie par 2 interVal pour avoir les 2 * 
00263             // les elem sym du bloc diagonal.
00264             fockFBB[lda][sma] += QC_TWO * interVal1;
00265             fockFBB[lda][sma] += interVal2;
00266           }
00267         }
00268       }
00269       // On remet la sous-matrice a sa place dans la matrice globale.
00270       this->matrix.insertSubMatrix(fockFAB,  firstAOofA, parameterA->getNbAO(),
00271                                    firstAOofB, parameterB->getNbAO() );
00272 
00273       // On remet la sous-matrice dans la matrice du systeme.
00274       this->matrix.insertSubTriangle(fockFBB, firstAOofB, parameterB->getNbAO() );
00275     }
00276     // On remet la sous-matrice dans la matrice du systeme.
00277     this->matrix.insertSubTriangle(fockFAA,  firstAOofA, parameterA->getNbAO() );
00278   }
00279   
00280   if (!manager.isDCComputation() && intgAcquisitionMethod == QC_INDIRECT_STORAGE &&  !isItfirstCall) {
00281     repInterElecIntegrals.getIntgReader()->rewindFile();
00282   }
00283 #ifdef VERBOSE_FOCK_BEFORE
00284   if (isItfirstCall) {
00285     ostringstream osstr;
00286     int id = workingSystem->getId();
00287     osstr << "fock_mine_" << id << "_before"; 
00288     this->matrix.printInFile(osstr.str().c_str());
00289   }
00290 #endif
00291 }
00295 template <class TPMatrix>
00296 template <class TPParam>
00297 void QCFockGtr<TPMatrix>::completeElems (const TPParam *    QCRestrict params,
00298                                          const QCSubDomain& QCRestrict workingDomain,
00299                                          const QCSubDomain& QCRestrict remoteDomain,
00300                                          const QCFloat *    QCRestrict remoteDensity,
00301                                          QCRepInterElec&    QCRestrict repInterElecIntegrals,
00302                                          QCSymMatrix&       QCRestrict fockFAA,
00303                                          const QCGeneralData&          data,
00304                                          int                           step,
00305                                          int                           nbDomains,
00306                                          QCIntgAcquisitionMethod       intgAcquisitionMethod,
00307                                          bool                          isFirstCall) 
00308 {
00309   
00310   // La matrice des integrales de repulsion
00311   QCMatrix& QCRestrict repIntegralsMatrix = repInterElecIntegrals.getMatrix();
00312 
00313   // les parametre de A et B
00314   const TPParam * QCRestrict paramA;
00315   const TPParam * QCRestrict paramB;
00316 
00317   // distance entre A et B
00318   QCFloat rAB;
00319 
00320   // Pour manipuler les coordonnees de A et B
00321   QCPoint3D coordsA;
00322   QCPoint3D coordsB;
00323 
00324   // La premiere orbitale de
00325   int firstAOofA;
00326 
00327   // Valeur intermediaire
00328   register QCFloat interVal1, interVal2;
00329   
00330   // les indices
00331   int i, j, k, mu, nu, lda, sma, muNu, ldaSma;
00332 
00333   // Les indices des block densite' 
00334   int blockIdx, firstBlockIdx;
00335 
00336   // la boucle sur les atomes du domaine courant
00337   for (i = 0; i < workingDomain.getNbAtoms(); ++i) {
00338     
00339     paramA     = &params[workingDomain.getParamIndexAt(i)];
00340     workingDomain.getPointAt(i, coordsA);
00341     firstAOofA = workingDomain.getFirstAOAt(i);
00342 
00343     this->matrix.extractSubTriangle(fockFAA, firstAOofA, paramA->getNbAO());
00344     
00345     // la boucle sur les atomes du domaine contribuant
00346     for (j = 0; j < remoteDomain.getNbAtoms(); ++j) {
00347       
00348       // On verifie si l'atome a deja apporte sa contribution
00349       bool contributed = false;
00350 
00351       if (data.getPartitionType() == QC_STANDARD_PART) {
00352         
00353         if (remoteDomain.getZone(j) == QC_SHELL2) {
00354           contributed = true;
00355           
00356         } else if (remoteDomain.getZone(j) == QC_SHELL1) {
00357           
00358           for (k=0; k <= step; ++k) {
00359             int contId    = (remoteDomain.getId() + k + 1) % nbDomains;
00360             int remoteIdx = remoteDomain.getRemoteIdx(contId, j);
00361             
00362             if (remoteIdx >= 0) {
00363               contributed = true;
00364               break;
00365             }
00366           }       
00367         }
00368                 
00369       } 
00370       else if ( data.getPartitionType() == QC_DIXON_PART ) {
00371         if (remoteDomain.getZone(j) != QC_CORE) {
00372           contributed = true;
00373         }
00374       }
00375       if (!contributed &&
00376           (workingDomain.getRemoteIdx(remoteDomain.getId(), i) != j)) {
00377   
00378         paramB     = &params[remoteDomain.getParamIndexAt(j)];
00379         remoteDomain.getPointAt(j, coordsB);
00380         //firstAOofB = remoteDomain.getFirstAOAt(j);
00381 
00382         rAB        = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00383                           QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00384                           QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00385 
00386 
00387         repInterElecIntegrals.computeElems(paramA,paramB, coordsA,coordsB, rAB,
00388                                            false, false,false, intgAcquisitionMethod,
00389                                            isFirstCall);
00390 
00391         firstBlockIdx = remoteDomain.getDensityBlockIdx(j);
00392         
00393         for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00394           muNu = mu * paramA->getNbAO();
00395           for (nu = 0;  nu <= mu; ++nu, ++muNu) {
00396             blockIdx   = firstBlockIdx;
00397             interVal1  = QC_ZERO;
00398             interVal2  = QC_ZERO;
00399             
00400             for (lda = 0; lda < paramB->getNbAO(); ++lda) {
00401               ldaSma = lda * paramB->getNbAO();
00402               for (sma = 0; sma < lda; ++sma, ++ldaSma, ++blockIdx) {
00403                 
00404                 // Completion des termes F_mu_mu ou F_mu_nu de l'atome A
00405                 // dans remoteDensity on a P_lda_sma
00406                 interVal1 += remoteDensity[blockIdx] * repIntegralsMatrix[muNu][ldaSma];
00407               }
00408 
00409               // Attention, ici, je profite du fait que sma == lda en sortie
00410               // de la bcle precedente, idem pour ldaSma.
00411               interVal2 +=  remoteDensity[blockIdx] * repIntegralsMatrix[muNu][ldaSma];
00412               ++blockIdx;
00413             }
00414             // On multiplie par 2 interVal pour avoir 2 * 
00415             // les elem sym du bloc diagonal.
00416             fockFAA[mu][nu] += QC_TWO * interVal1;
00417             fockFAA[mu][nu] += interVal2;
00418           }
00419         }
00420       }          
00421     }
00422 
00423     // On remet la sous-matrice dans la matrice du systeme.
00424     this->matrix.insertSubTriangle(fockFAA, firstAOofA, paramA->getNbAO());
00425   }
00426 }
00431 template <class TPMatrix>
00432 template <class TPIterator>
00433 void
00434 QCFockGtr<TPMatrix>::levelShifting (TPIterator& workingSystem, const QCFloat& doubleGapShift) {
00435   //
00436   // La matrice densite
00437   QCSymMatrix& QCRestrict densityMatrix = workingSystem->getMatrices().getDensityP().getMatrix();
00438 
00439   // Le shifter
00440   QCFloat gapShift = QC_HALF * doubleGapShift;
00441   int i, j;
00442   //
00443   // On commence par le debut, soit la completions des termes diagonaux F_mu_mu
00444   //
00445   for (i = 0; i < densityMatrix.getDim() ; ++i) {
00446     for (j = 0; j <= i; ++j) {
00447       // gapShift est negatif, on envoie les orb occ vers le bas.
00448       this->matrix[i][j] += gapShift * densityMatrix[i][j];
00449     }
00450     // Puis on remonte tout le monde
00451     this->matrix[i][i] -= doubleGapShift;
00452     // Au final, on a ecarte les virtuelles des occupees.
00453     //   Les virtuelles ont ete remontees de doubleGapShift.
00454     // C'est la valeur qui manquait au gap pour atteindre la valeur
00455     // de gap minimum definie par le shiftingLevelParameter
00456   }
00457 }
00458 
00465 template <class TPMatrix>
00466 template <class TPManager>
00467 void
00468 QCFockGtr<TPMatrix>::diagonalize (TPManager& manager,
00469                                   typename TPManager::TSystem::QCIterator& workingSystem,
00470                                   const int& nbIter) {
00471 
00472   // les gestionnaire de memoire
00473   QCMemory& memory = manager.getMemory();
00474 
00475   // Les matrices courantes
00476   QCModelMatrices& QCRestrict matrices = 
00477     workingSystem->getMatrices();
00478 
00479   // La matrice des vecteurs propres
00480   QCMatrix& transEigenVectMatrix = matrices.getEigenVectCt().getMatrix();
00481 
00482   // Les valeures propres
00483   QCFloat * QCRestrict eigenVal = matrices.getEigenVal();
00484 
00485   // Le tableau des degenerescence.
00486   int * QCRestrict degener         = matrices.getDegener();
00487   const QCFloat    toleranceFactor = 1e-12;
00488 
00489   // L'algo de diagonalisation
00490   QCDiagoAlgorithm diagoAlgorithm = 
00491     manager.getGeneralData().getDiagoAlgorithm();
00492 
00493 #ifdef USE_LAPACK
00494   // Si on utilise une diago Divide and Conquer, la matrice des vecteurs propres
00495   // contient en entree de la diago, la matrice a diagonaliser.
00496   if (diagoAlgorithm  == QC_DC_DIAGO) {
00497 
00498     // On copie la matrice de Fock dans celle des vecteurs propres 
00499     // qui sera a la fois l'entree et la sortie de la subroutine de 
00500     // diagonalisation fournie par LAPACK.
00501     // La methode copyIn a ete situee dans QCSymMatrix et non dans QCMatrix
00502     // car QCMatrix.h est inclus dans QCSymMatrix.hpp et le contraire est
00503     // impossible.
00504     this->matrix.copyIn(transEigenVectMatrix);
00505   }
00506 #endif // USE_LAPACK
00507   
00508   // La diagonalisation.
00509   if (!this->matrix.diagonalize(memory, diagoAlgorithm, toleranceFactor, eigenVal,
00510                           degener,transEigenVectMatrix.getElems(), nbIter)) {
00511     cerr << __FILE__ << "diagonalise: diagonalisation aborted"
00512          << endl;
00513     exit(EXIT_FAILURE);
00514   }
00515 }
00516 
00517 
00518 
00519 
00520 
00521 
00522 
00526 template class QCFockGtr<QCSymMatrix>;
00527 
00528 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::completeElems,
00529                                     const bool&);
00530 
00531 QCPARAMETER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::completeElems,
00532                                  TEN_PARAMS(const QCSubDomain& QCRestrict,
00533                                             const QCSubDomain& QCRestrict,
00534                                             const QCFloat *    QCRestrict,
00535                                             QCRepInterElec&    QCRestrict,
00536                                             QCSymMatrix&       QCRestrict,
00537                                             const QCGeneralData&,
00538                                             int,
00539                                             int,
00540                                             QCIntgAcquisitionMethod,
00541                                             bool));
00542 
00543 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::diagonalize,
00544                                     const int&);
00545 
00546 QCITERATOR_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::levelShifting,
00547                                 const QCFloat&);
00548 

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