QCHamiltonGtr.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 
00021 
00022 #include "QCCommon.hpp"
00023 #include "QCHamiltonGtr.hpp"
00024 #include "QCSymMatrix.hpp"
00025 #include "QCRepInterElec.hpp"
00026 #include "QCManager.hpp"
00027 #include "QCMndo.hpp"
00028 #include "QCAm1.hpp"
00029 #include "QCPm3.hpp"
00030 #include "QCMndoParam.hpp"
00031 #include "QCAm1Param.hpp"
00032 #include "QCPm3Param.hpp"
00033 #include "QCGlobalSystem.hpp"
00034 #include "QCMDSystem.hpp"
00035 #include "QCDistMDSystem.hpp"
00036 #include "QCSCFAlgo.hpp"
00037 #include "QCDCAlgo.hpp"
00038 
00039 
00040 
00044 template <class TPMatrix>
00045 QCHamiltonGtr<TPMatrix>::QCHamiltonGtr (void) : 
00046   QCMatElemGtr<TPMatrix>() 
00047 {}
00048 
00049 
00050 
00054 template <class TPMatrix>
00055 QCHamiltonGtr<TPMatrix>::~QCHamiltonGtr (void) 
00056 {}
00057 
00058 
00059 
00060 
00064 template <class TPMatrix>
00065 template <class TPManager>
00066 void
00067 QCHamiltonGtr<TPMatrix>::computeElems (TPManager& manager,
00068                                        typename TPManager::TSystem::QCIterator& workingSystem) {
00069 
00070   
00071   typedef typename TPManager::TModel         TModel;
00072   typedef typename TPManager::TModel::TParam TParam;
00073   
00078   const QCSystem& QCRestrict currentSystem = workingSystem->getSystem();
00079   TModel&         QCRestrict model         = manager.getModel();
00080   const TParam *  QCRestrict params        = manager.getParameters();
00081 
00082   //
00083   const QCGeneralData& data = manager.getGeneralData();
00084 
00085   // les parametres de A et de B
00086   const TParam * QCRestrict paramA;
00087   const TParam * QCRestrict paramB;
00088 
00089   // distance entre A et B
00090   QCFloat rAB;
00091 
00092   // La matrice d'overlap
00093   QCSymMatrix& QCRestrict overlapMatrix = 
00094     workingSystem->getMatrices().getOverlapS()->getMatrix();
00095 
00096   // Le generateur de repultion inter elec
00097   QCRepInterElec& QCRestrict repInterElecIntegrals =
00098     workingSystem->getMatrices().getRepInterElecIntegrals();
00099 
00100   // La matrice des integrales de repulsion
00101   QCMatrix& QCRestrict repIntegralsMatrix =
00102     repInterElecIntegrals.getMatrix();
00103 
00104   // matrices de Hamilton locales a cette methode. Attention,
00105   // cette matrice n'est pas symetrique
00106   QCSymMatrix& QCRestrict hamiltonHAA = model.getSpWorkingAA();
00107   QCMatrix&    QCRestrict hamiltonHAB = model.getSpWorkingAB();
00108 
00109   // matrice Sab locale a cette methode.
00110   QCMatrix&    QCRestrict overlapSAB  = model.getSpdWorkingAB();
00111 
00112   // l'index de ss dans le tableau des integrales de rep inter elec.
00113   const int  SSRepIndex = 0;
00114 
00115   // Les premiers index d'orbitales de A et B dans la matrice complete du
00116   // systeme et l'index du premier index du successeur de A tjs dans la 
00117   // matrice complete.
00118   int firstAOofA, firstAOofB;
00119 
00120   // Pour manipuler les coordonnees
00121   QCPoint3D coordsA, coordsB;
00122 
00123   // les indices
00124   int i, j, mu, nu, lda, muNu;
00125   
00126   // la double boucle sur les atomes
00127   for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00128     
00129     paramA     = &params[currentSystem.getParamIndexAt(i)];
00130     firstAOofA = currentSystem.getFirstAOAt(i);
00131     currentSystem.getPointAt(i, coordsA);
00132 
00133     // Je mets les elements a 0 car on fait des += a un moment.
00134     hamiltonHAA.initClean();
00135 
00136     // On remplit les termes diagonaux : H_mu_mu, 
00137     // on ajoute Uss et Usp.
00138     for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00139 
00140       // Attention, en l'occurence, il s agit de hamiltonHAA.
00141       hamiltonHAA[mu][mu] = paramA->getU(mu);
00142     }
00143     
00144     // Boucle sur les B du meme systeme pour remplir les H_mu_nu.
00145     for (j = 0; j < currentSystem.getNbAtoms(); ++j) {
00146       if (j != i) {
00147 
00148         if ( (data.getPartitionType() == QC_STANDARD_PART && 
00149               workingSystem->getZone(j) != QC_SHELL2) || 
00150              
00151              (data.getPartitionType() == QC_DIXON_PART &&
00152               workingSystem->getZone(j) == QC_CORE) ) {
00153           
00154 //        if (workingSystem->getId() == 0 && i == 14) {
00155 //          cout << i 
00156 //               << " <- " << j << " (" << workingSystem->getGlobalIdx(j) << ")" 
00157 //               << endl; 
00158 //        }
00159           
00160 
00161           paramB  = &params[currentSystem.getParamIndexAt(j)];
00162           currentSystem.getPointAt(j, coordsB);
00163           rAB     = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00164                          QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00165                          QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00166           
00167           repInterElecIntegrals.computeElems(paramB,
00168                                              paramA,
00169                                              coordsB,
00170                                              coordsA,
00171                                              rAB,
00172                                              true,
00173                                              false);
00174           
00175           
00176           // On remplit les termes H_mu_nu.
00177           // k prend le nb d'OA de A valeurs au carre.
00178           for (mu = 0, muNu = 0; mu < paramA->getNbAO(); ++mu) {
00179             muNu = mu * paramA->getNbAO();
00180             for (nu = 0; nu <= mu; ++nu, ++muNu) {
00181               hamiltonHAA[mu][nu] -= 
00182                 paramB->getCoreCharge() * repIntegralsMatrix[SSRepIndex][muNu];
00183             }
00184           }
00185         }
00186       } 
00187     }      
00188 
00189     
00190     
00191     // et on insere la sous matrice dans la matrice principale.
00192     // Attention, il s agit de la sous-matrice hamiltonHAA en fait.
00193     this->matrix.insertSubTriangle(hamiltonHAA,
00194                              firstAOofA, paramA->getNbAO());
00195 
00196     
00197     // Le cas ou A est different de B.
00198     // On veut assigner les H_mu_lambda.
00199     for (j = 0; j < i; ++j) {
00200       
00201       paramB     = &params[currentSystem.getParamIndexAt(j)];
00202       firstAOofB = currentSystem.getFirstAOAt(j);
00203       
00204       // On extrait la sous matrice d'overlap dont on a besoin.
00205       overlapMatrix.extractSubMatrix(overlapSAB,
00206                                      firstAOofA, paramA->getNbAO(),
00207                                      firstAOofB, paramB->getNbAO());
00208       
00209       // remplissage de la matrice de Hamilton
00210       // avec f4 * S_mu_lambda.
00211       for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00212         for (lda = 0; lda < paramB->getNbAO(); ++lda) {
00213 
00214           // Ok, c est une assignation pure, on ne remet pas a 0.
00215           hamiltonHAB[mu][lda] = 
00216             ( (paramA->getBeta(mu) + paramB->getBeta(lda) ) * 
00217               QC_HALF * overlapSAB[mu][lda] );
00218         }
00219       }
00220       
00221       // On insere la sous-matrice a sa place dans la matrice globale.
00222       this->matrix.insertSubMatrix(hamiltonHAB, 
00223                              firstAOofA, paramA->getNbAO(),
00224                              firstAOofB, paramB->getNbAO() );
00225       
00226     }
00227   }
00228 #ifdef VERBOSE_HAMILTON_BEFORE
00229   ostringstream osstr;
00230   int id = workingSystem->getId();
00231   osstr << "hamilton_mine_";
00232   if (id >= 0) {
00233     osstr << id << "_before";
00234   }
00235   this->matrix.printInFile(osstr.str().c_str());
00236 #endif
00237 }
00238 
00239 
00240 
00241 
00242 
00246 template <class TPMatrix>
00247 template <class TPParam>
00248 void
00249 QCHamiltonGtr<TPMatrix>::completeElems (const TPParam *    QCRestrict params,
00250                                         const QCSubDomain& QCRestrict workingDomain,
00251                                         const QCSubDomain& QCRestrict remoteDomain,
00252                                         QCRepInterElec&    QCRestrict repInterElecIntegrals,
00253                                         QCSymMatrix&       QCRestrict hamiltonHAA,
00254                                         const QCGeneralData&          data,
00255                                         int                           step,
00256                                         int                           nbDomains) 
00257 {
00258   
00259 
00260   // La matrice des integrales de repulsion
00261   QCMatrix& QCRestrict repIntegralsMatrix =
00262     repInterElecIntegrals.getMatrix();
00263 
00264   // les parametre de A et B
00265   const TPParam * QCRestrict paramA;
00266   const TPParam * QCRestrict paramB;
00267 
00268   // distance entre A et B
00269   QCFloat rAB;
00270 
00271   // l'index de ss dans le tableau des integrales de rep inter elec.
00272   const int  SSRepIndex = 0;
00273 
00274   // Pour manipuler les coordonnees de A et B
00275   QCPoint3D coordsA;
00276   QCPoint3D coordsB;
00277 
00278   // les indices
00279   int i, j, k, mu, nu, muNu;
00280 
00281   // la boucle sur les atomes du domaine courant
00282   for (i = 0; i < workingDomain.getNbAtoms(); ++i) {
00283 
00284     paramA  = &params[workingDomain.getParamIndexAt(i)];
00285     workingDomain.getPointAt(i, coordsA);
00286 
00287 //     matrix.extractSubTriangle(hamiltonHAA,
00288 //                            workingDomain.getFirstAOAt(i), 
00289 //                            paramA->getNbAO());
00290 
00291     hamiltonHAA.initClean();
00292 
00293     // la boucle sur les atomes du domaine contribuant
00294     for (j = 0; j < remoteDomain.getNbAtoms(); ++j) {
00295       
00296       // On verifie si l'atome a deja apporte sa contribution
00297       bool contributed = false;
00298       
00299       if (data.getPartitionType() == QC_STANDARD_PART) {
00300         
00301         if (remoteDomain.getZone(j) == QC_SHELL2) {
00302           contributed = true;
00303           
00304         } else if (remoteDomain.getZone(j) == QC_SHELL1) {
00305           
00306           for (k=0; k <= step; ++k) {
00307             int contId    = (remoteDomain.getId() + k + 1) % nbDomains;
00308             int remoteIdx = remoteDomain.getRemoteIdx(contId, j);
00309             
00310             if (remoteIdx >= 0) {
00311               contributed = true;
00312               break;
00313             }
00314           }       
00315         }
00316                 
00317       } else if ( data.getPartitionType() == QC_DIXON_PART ) {
00318         
00319         if (remoteDomain.getZone(j) != QC_CORE) {
00320           contributed = true;
00321         }
00322       }
00323       
00324       
00325       if (!contributed &&
00326           (workingDomain.getRemoteIdx(remoteDomain.getId(), i) != j)) {
00327 
00328 //      if (workingDomain.getId() == 0 &&i == 14) {
00329 //        cout << i 
00330 //             << " <- " << j << " (" << remoteDomain.getGlobalIdx(j) << ")" 
00331 //             << endl; 
00332 //      }
00333 
00334         paramB  = &params[remoteDomain.getParamIndexAt(j)];
00335         remoteDomain.getPointAt(j, coordsB);
00336         rAB     = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00337                        QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00338                        QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00339         
00340         repInterElecIntegrals.computeElems(paramB,
00341                                            paramA,
00342                                            coordsB,
00343                                            coordsA,
00344                                            rAB,
00345                                            true,
00346                                            false);
00347 
00348         // On remplit les termes H_mu_nu.
00349         // k prend le nb d'OA de A valeurs au carre.
00350         for (mu = 0, muNu = 0; mu < paramA->getNbAO(); ++mu) {
00351           muNu = mu * paramA->getNbAO();
00352           for (nu = 0; nu <= mu; ++nu, ++muNu) {
00353             hamiltonHAA[mu][nu] -= 
00354               paramB->getCoreCharge() * repIntegralsMatrix[SSRepIndex][muNu];
00355           }
00356         }
00357       }
00358     }
00359 
00360     // et on insere la sous matrice dans la matrice principale.
00361     // Attention, il s agit de la sous-matrice hamiltonHAA en fait.
00362     this->matrix.completeSubTriangle(hamiltonHAA,
00363                                workingDomain.getFirstAOAt(i), 
00364                                paramA->getNbAO());
00365   }
00366 
00367 }
00368 
00369 
00370 
00371 
00372 
00373 
00374 
00375 
00379 template class QCHamiltonGtr<QCSymMatrix>;
00380 QCMANAGER_ITER_METH_EXPL_INST(void QCHamiltonGtr<QCSymMatrix>::computeElems);
00381 
00382 QCPARAMETER_METH_EXPL_INST_PARAM(void QCHamiltonGtr<QCSymMatrix>::completeElems,
00383                                  SEVEN_PARAMS(const QCSubDomain& QCRestrict,
00384                                               const QCSubDomain& QCRestrict,
00385                                               QCRepInterElec& QCRestrict,
00386                                               QCSymMatrix&    QCRestrict,
00387                                               const QCGeneralData&,
00388                                               int,
00389                                               int));

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