QCOverlapGtr.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 "QCCommon.hpp"
00021 #include "QCOverlapGtr.hpp"
00022 #include "QCSymMatrix.hpp"
00023 #include "QCManager.hpp"
00024 #include "QCMndo.hpp"
00025 #include "QCMndoParam.hpp"
00026 #include "QCAm1.hpp"
00027 #include "QCAm1Param.hpp"
00028 #include "QCPm3.hpp"
00029 #include "QCPm3Param.hpp"
00030 #include "QCGlobalSystem.hpp"
00031 #include "QCMDSystem.hpp"
00032 #include "QCDistMDSystem.hpp"
00033 #include "QCSCFAlgo.hpp"
00034 #include "QCDCAlgo.hpp"
00035 
00036 
00037 
00038 //#define OUTPUT_DERIVATIVES
00039 //#define OUTPUT_DERIVATIVES2
00040 
00041 
00042 
00046 template <class TPMatrix>
00047 QCFloat *QCOverlapGtr<TPMatrix>::ExpHalfVector = NULL;
00048 
00049 
00050 
00054 template <class TPMatrix>
00055 QCOverlapGtr<TPMatrix>::QCOverlapGtr (void) : 
00056   QCMatElemGtr<TPMatrix> () 
00057 {}
00058 
00059 
00060 
00061 
00065 template <class TPMatrix>
00066 QCOverlapGtr<TPMatrix>::QCOverlapGtr (const int& dim) : 
00067   QCMatElemGtr<TPMatrix> (dim) 
00068 {}
00069 
00070 
00071 
00072 
00076 template <class TPMatrix>
00077 QCOverlapGtr<TPMatrix>::~QCOverlapGtr (void) 
00078 {}
00079 
00080 
00081 
00082 
00086 template <class TPMatrix>
00087 template <class TPManager>
00088 void    
00089 QCOverlapGtr<TPMatrix>::FillExpHalfVector (TPManager& manager) 
00090 {
00091 
00092   static int*     quantNumNVector;
00093   static QCFloat* dzetaVector;
00094   
00095   typedef typename TPManager::TModel TModel;
00096 
00097   const typename TModel::TParam *QCRestrict params = manager.getParameters();
00098 
00099   // Des tableaux utilistaires pour l appel a overlp.
00100   quantNumNVector = new int     [TModel::NbMaxAtomTypes];
00101   dzetaVector     = new QCFloat [TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes];
00102   ExpHalfVector   = new QCFloat [TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes];
00103 
00104   // On n oublie pas d'initialiser les tableaux avec des val nulles.
00105   memset(quantNumNVector, 0, 
00106          TModel::NbMaxAtomTypes * sizeof(int) );
00107   memset(dzetaVector,     0,
00108          TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes * sizeof(QCFloat) );
00109   memset(ExpHalfVector,   0,
00110          TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes * sizeof(QCFloat) );
00111 
00112   // On remplit avec les nb quant princ et les dzeta trouves dans les classes 
00113   // de param.
00114   FillQuantNumNVector(params, manager.getNbAtomTypes(), quantNumNVector);
00115   FillDzetaVector    (params, manager.getNbAtomTypes(), 
00116                       TModel::NbMaxOrbTypes, dzetaVector);  
00117 
00123   QCFortran(fillexhalf)
00124     (quantNumNVector, dzetaVector, ExpHalfVector);
00125   
00126   if (quantNumNVector) { delete [] quantNumNVector; }
00127   if (dzetaVector)     { delete [] dzetaVector; }
00128 }
00129 
00130 
00131 
00132 
00136 template <class TPMatrix>
00137 void    
00138 QCOverlapGtr<TPMatrix>::DestroyExpHalfVector (void) {
00139   if (ExpHalfVector) 
00140     { 
00141       delete [] ExpHalfVector; 
00142       ExpHalfVector = NULL; 
00143     }
00144 }
00145 
00146 
00147 
00148 
00152 template <class TPMatrix>
00153 template <class TPManager>
00154 void    
00155 QCOverlapGtr<TPMatrix>::computeElems (TPManager& manager,
00156                                       typename TPManager::TSystem::QCIterator& workingSystem,
00157                                       // Le parametre de shift utile si
00158                                       // on fait des derivations
00159                                       // numeriques.
00160                                       const QCPoint3D& shift) 
00161 {
00162   
00163   typedef typename TPManager::TModel         TModel;
00164   typedef typename TPManager::TModel::TParam TParam;
00165 
00170   TModel&        QCRestrict model  = manager.getModel();
00171   const TParam * QCRestrict params = manager.getParameters(); 
00172 
00177   QCMatrix& QCRestrict overlapSAB = model.getSpdWorkingAB();
00178 
00179 
00180   computeElems(params, 
00181                workingSystem->getSystem(),
00182                overlapSAB,
00183                TModel::NbMaxOrbTypes,
00184                shift);
00185 }
00186 
00187 
00188 
00189 
00190 
00194 template <class TPMatrix>
00195 template <class TPParam>
00196 void 
00197 QCOverlapGtr<TPMatrix>::computeElems (const TPParam *  QCRestrict parameters,
00198                                       const QCSystem&  QCRestrict system,
00199                                       QCMatrix&        QCRestrict overlapSAB,
00200                                       int              nbMaxOrbTypes,
00201                                       const QCPoint3D& QCRestrict shift) {
00202   
00207   int   firstAOofA, firstAOofB;
00208 
00212   const TPParam *QCRestrict paramA;
00213   const TPParam *QCRestrict paramB;
00214 
00215   
00216   QCPoint3D coordsA, coordsB;
00217 
00218 
00223   for (int i = 1; i < system.getNbAtoms(); ++i) {
00224 
00225     system.getPointAt(i, coordsA);
00226     paramA = &parameters[system.getParamIndexAt(i)];
00227     firstAOofA = system.getFirstAOAt(i);
00228     
00229     for (int mu = firstAOofA; 
00230          mu < (firstAOofA + paramA->getNbAO()); 
00231          mu++) 
00232       {
00233         // Je mets des 1 sur la diagonale.
00234         this->matrix[mu][mu] = QC_ONE;
00235       }
00236 
00237     for (int j = 0; j < i; j++) {
00238 
00239       system.getPointAt(j, coordsB);
00240       paramB = &parameters[system.getParamIndexAt(j)];
00241       firstAOofB = system.getFirstAOAt(j);
00242       
00243       computeElemsAux ( coordsA,
00244                         *paramA,
00245                         (paramA->getAtomicNbZ()-1) * nbMaxOrbTypes,
00246                         coordsB,
00247                         *paramB,
00248                         (paramB->getAtomicNbZ()-1) * nbMaxOrbTypes,
00249                         overlapSAB,
00250                         shift );
00251       
00258       this->matrix.insertSubMatrix( overlapSAB, 
00259                               firstAOofA, paramA->getNbAO(),
00260                               firstAOofB, paramB->getNbAO() );
00261       
00262     }
00263   }
00264 #ifdef DEBUG
00265   //matrix.print (">>>>>>>>>> OVERLAP MATRIX <<<<<<<<<<<<<");
00266 #endif
00267 }
00268 
00269 
00270 
00271 
00272 
00276 template <class TPMatrix>
00277 template <class TPParam>
00278 void
00279 QCOverlapGtr<TPMatrix>::computeElemsAux (const QCPoint3D& QCRestrict coordsA,
00280                                          const TPParam&   QCRestrict paramA,
00281                                          int              expHalfIndexA,
00282                                          const QCPoint3D& QCRestrict coordsB,
00283                                          const TPParam&   QCRestrict paramB,
00284                                          int              expHalfIndexB,
00285                                          QCMatrix&        QCRestrict overlapSAB,
00286                                          const QCPoint3D& shift) 
00287 {
00288   
00289   
00290   QCPoint3D coordsdB;
00291 
00292   coordsdB[COORDX] = coordsB[COORDX] + shift[COORDX];
00293   coordsdB[COORDY] = coordsB[COORDY] + shift[COORDY];
00294   coordsdB[COORDZ] = coordsB[COORDZ] + shift[COORDZ];
00295 
00296   QCFloat rAB = sqrt( QCPow<2>(coordsA[COORDX] - coordsdB[COORDX]) +
00297                       QCPow<2>(coordsA[COORDY] - coordsdB[COORDY]) +
00298                       QCPow<2>(coordsA[COORDZ] - coordsdB[COORDZ]) );
00299   
00300   QCFortran(overlp)
00301     (//
00302      // Attention, l indice du premier entree est celui
00303      // qui varie le plus vite.
00304      paramB.getAtomicNbZ(),
00305      paramB.getNbAO(),
00306      paramB.getQuantumNumberN(),
00307      paramB.getDzetaS(),
00308      paramB.getDzetaP(),
00309      paramB.getDzetaD(),
00310      coordsdB[COORDX],
00311      coordsdB[COORDY],
00312      coordsdB[COORDZ],
00313      &ExpHalfVector[expHalfIndexB],
00314      paramA.getAtomicNbZ(),
00315      paramA.getNbAO(),
00316      paramA.getQuantumNumberN(),
00317      paramA.getDzetaS(),
00318      paramA.getDzetaP(),
00319      paramA.getDzetaD(),
00320      coordsA[COORDX],
00321      coordsA[COORDY],
00322      coordsA[COORDZ],
00323      &ExpHalfVector[expHalfIndexA],
00324      rAB,
00325      overlapSAB.getElems() 
00326      );
00327   
00328 }
00329 
00330 
00331 
00332 
00333 
00337 template <class TPMatrix>
00338 template <class TPManager>
00339 void  
00340 QCOverlapGtr<TPMatrix>::computeDerivatives (TPManager&       manager,
00341                                             typename TPManager::TSystem::QCIterator& workingSystem,
00342                                             const QCPoint3D& shift,
00343                                             const QCFloat&   QCRestrict inv2DeltaQ) {
00344 
00345   QCOverlapGtr<QCSymMatrix>& QCRestrict 
00346     workingIntg = *(workingSystem->getMatrices().getOverlapS());
00347   
00348 
00349   if (manager.getGeneralData().getDerivationType() == QC_NUMERICAL_DERIVATION) {
00350     
00351     //On fait le calcul des integrales de rep inter elec
00352     //pour r-deltaQ
00353     
00354     QCPoint3D mshift;
00355     mshift -= shift;
00356     workingIntg.computeElems(manager, workingSystem, mshift);
00357 
00358 #ifdef OUTPUT_DERIVATIVES2
00359     workingIntg.getMatrix().printInFile("DERIVATIVES_MS");
00360 #endif
00361 
00362     
00363     //On fait le calcul des integrales de rep inter elec
00364     //pour r+deltaQ
00365     computeElems(manager, workingSystem, shift);
00366 
00367 #ifdef OUTPUT_DERIVATIVES2
00368     this->matrix.printInFile("DERIVATIVES_S");
00369 #endif
00370     
00371     //Et pour finir, on calcule les derivees.
00372     this->matrix.calculateNumericalDerivatives(workingIntg.getMatrix(),
00373                                          inv2DeltaQ);
00374    
00375    
00376   } else {//manager.getGeneralData().getDerivationType() == QC_ANALYTICAL_DERIVATION 
00377     
00378   }
00379 
00380 #ifdef OUTPUT_DERIVATIVES
00381     matrix.printInFile("DERIVATIVES");
00382 #endif
00383 
00384 }
00385 
00386 
00387 
00388 
00389 
00390 
00391 /*************************
00392  *       PRIVATE         *
00393  *************************/
00394 
00395 
00401 template <class TPMatrix>
00402 template <class TPParam>
00403 void
00404 QCOverlapGtr<TPMatrix>::FillQuantNumNVector (const TPParam *QCRestrict params,
00405                                              const int& nbParameters,
00406                                              int* QCRestrict quantNumNVector) 
00407 {
00408   for (int i=0; i < nbParameters; ++i) {
00409     quantNumNVector[params[i].getAtomicNbZ()-1] = params[i].getQuantumNumberN();
00410   }
00411 }
00412 
00413 
00414 
00415 
00421 template <class TPMatrix>
00422 template <class TPParam>
00423 void
00424 QCOverlapGtr<TPMatrix>::FillDzetaVector (const TPParam *QCRestrict params,
00425                                          const int& nbParameters,
00426                                          const int& nbMaxOrbTypes,
00427                                          QCFloat* QCRestrict dzetaVector) 
00428 {
00429   int index;
00433   for (int i=0; i < nbParameters; ++i) {
00434     
00435     //On calcule l index a partir duquel vont se trouver les 3 dzeta.
00436     index = (params[i].getAtomicNbZ()-1) * nbMaxOrbTypes;
00437     
00438     //en fait les dzetaS, dzetaP et dzetaD de l atome sont 
00439     //stocke continuement dans le tableau.
00440     dzetaVector[index] = params[i].getDzetaS();
00441     dzetaVector[index+1] = params[i].getDzetaP();
00442     
00443     //Attention, la fonction getDzetaD retourne 0e0 pour le moment.
00444     dzetaVector[index+2] = params[i].getDzetaD();
00445   }
00446 }
00447 
00448 
00449 
00450 
00454 template class QCOverlapGtr<QCSymMatrix>;
00455 
00456 QCMANAGER_METH_EXPL_INST(void QCOverlapGtr<QCSymMatrix>::FillExpHalfVector);
00457 
00458 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCOverlapGtr<QCSymMatrix>::computeElems,
00459                                     const QCPoint3D&);
00460 
00461 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCOverlapGtr<QCSymMatrix>::computeDerivatives,
00462                                     TWO_PARAMS(const QCPoint3D&,
00463                                                const QCFloat& QCRestrict));

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