QCMndo.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 "QCMndo.hpp"
00021 #include "QCAm1.hpp"
00022 #include "QCPm3.hpp"
00023 #include "QCMndoParam.hpp"
00024 #include "QCAm1Param.hpp"
00025 #include "QCPm3Param.hpp"
00026 #include "QCGlobalSystem.hpp"
00027 #include "QCMDSystem.hpp"
00028 #include "QCDistMDSystem.hpp"
00029 #include "QCSCFAlgo.hpp"
00030 #include "QCDCAlgo.hpp"
00031 #include "QCManager.hpp"
00032 #include "QCIntgReader.hpp"
00033 #include "QCIntgWriter.hpp"
00034 
00035 
00036 
00037 #ifndef QC_NO_DEBUG
00038 #  define VERBOSE_DERIV2
00039 #endif
00040 
00041 
00042 
00047 const int   QCMndo::NbMaxOrbTypes = 3;
00048 
00049 
00054 const int   QCMndo::NbMaxAtomTypes = 83;
00055 
00056 
00057 
00061 QCMndo::QCMndo (const int& nbAOsp, const int& nbAOspd) :
00062   QCModel(nbAOsp, nbAOspd)
00063 {}
00064 
00065 
00066 
00067 
00071 QCMndo::~QCMndo (void)
00072 {
00073 
00074 }
00075 
00076 
00077 
00078 
00082 template <class TPManager>
00083 QCFloat
00084 QCMndo::computeCoreEnergy (TPManager& manager,
00085                            typename TPManager::TSystem::QCIterator& workingSystem) {
00086 
00087   // On recupere avec des restricts pour le calcul les
00088   // structures qui contiennent les donnees a manipuler
00089   QCSystem& QCRestrict currentSystem = workingSystem->getSystem();
00090   QCFloat coreEnergy;
00091   
00092   // Le modele
00093   typename TPManager::TModel& QCRestrict model = manager.getModel();
00094 
00095   // Les matrices
00096   QCModelMatrices& QCRestrict matrices = workingSystem->getMatrices();
00097   QCRepInterElec&  QCRestrict repInterElecIntegrals = matrices.getRepInterElecIntegrals();
00098   
00099   // l'index de ss dans le tableau des integrales de rep inter elec.
00100   const int SSRepIndex = 0;
00101 
00102   // Les parametres
00103   const typename TPManager::TParam *QCRestrict parameters =
00104     manager.getParameters();
00105   const typename TPManager::TParam *QCRestrict parameterA;
00106   const typename TPManager::TParam *QCRestrict parameterB;
00107 
00108   // Les coords
00109   QCPoint3D coordsA;
00110   QCPoint3D coordsB;
00111 
00112   // Les types d atomes
00113   int atomicNbA;
00114   int atomicNbB;
00115 
00116 
00117   // Distance entre a et b.
00118   QCFloat   rAB;  
00119   QCFloat   expMalfAR, expMalfBR;
00120   int i, j;
00121 
00122   // La double boucle sur les atomes.
00123   // je prends i < j.
00124   coreEnergy = QC_ZERO;
00125   for (i = 1; i < currentSystem.getNbAtoms(); ++i) {
00126 
00127     currentSystem.getPointAt(i, coordsA);
00128     parameterA = &parameters[currentSystem.getParamIndexAt(i)];
00129     atomicNbA  = parameterA->getAtomicNbZ();
00130     
00131     for (j = 0; j < i; ++j) {
00132 
00133       currentSystem.getPointAt(j, coordsB);
00134       parameterB = &parameters[currentSystem.getParamIndexAt(j)];
00135       atomicNbB  = parameterB->getAtomicNbZ();
00136       rAB        = currentSystem.distance(i,j);
00137       expMalfAR  = QCExp(-parameterA->getAlpha() * rAB);
00138       expMalfBR  = QCExp(-parameterB->getAlpha() * rAB);
00139 
00140       // Calcul des elements de la matrice de rep inter elec
00141       repInterElecIntegrals.computeElems(parameterA,
00142                                          parameterB,
00143                                          coordsA,
00144                                          coordsB,
00145                                          rAB,
00146                                          //ssACalc recoit true,
00147                                          true,
00148                                          //
00149                                          //ssBCalc recoit true,
00150                                          //car ici, on ne fait les calculs que pour 
00151                                          //l'element <ss|ss> de la matrice de rep
00152                                          //inter-elec.
00153                                          true);
00154 
00155       coreEnergy += 
00156         parameterA->getCoreCharge() * 
00157         parameterB->getCoreCharge() *
00158         (repInterElecIntegrals.getMatrix()[SSRepIndex][SSRepIndex] *
00159          (QC_ONE +
00160           model.sometimesR(atomicNbA, atomicNbB, rAB) * expMalfAR +
00161           model.sometimesR(atomicNbB, atomicNbA, rAB) * expMalfBR) +
00162          model.exponentialExpansion(parameterA, parameterB, rAB) / rAB);
00163     }
00164   }
00165   
00166   workingSystem->setCoreEnergy(coreEnergy);
00167 
00168   return coreEnergy;
00169 }
00170 
00171 
00172 
00173 
00177 template <class TPIterator>
00178 QCFloat
00179 QCMndo::computeElecEnergy (TPIterator& workingSystem,
00180                            QCDensityGtr<QCSymMatrix>& QCRestrict densityP) {
00181   QC_TRACE_ENER(" BEGIN QCMndo::computeElecEnergy ");
00182   // L'energie
00183   QCFloat elecEnergy = QC_ZERO;
00184 
00185   // L'energie electronique est egale a 1/2 tr(P(H+F) ).
00186   elecEnergy = ( QC_HALF * 
00187                  ( workingSystem->getHamiltonH().traceProduct (densityP) + 
00188                    workingSystem->getFockF().traceProduct     (densityP) ) );
00189   
00190 if (finite(elecEnergy)) {
00191 
00192     workingSystem->setElecEnergy(elecEnergy);
00193 
00194   } else {
00195 
00196 #ifdef DEV_SCF_CLASS_NO_MEMORY
00197 
00203     densityP.swapMatrix(workingSystem->getDensityPnm1());
00204 #endif // DEV_SCF_CLASS_NO_MEMORY
00205     
00209     cerr << __FILE__ << " computeElecEnergy: bad value of electronic energy : "<<elecEnergy
00210          << endl;
00211     exit(EXIT_FAILURE);
00212   }
00213  QC_TRACE_ENER("END   QCMndo::computeElecEnergy "<<elecEnergy);
00214 
00215   return elecEnergy;
00216 }
00217 
00218 
00219 
00220 
00221 
00222 
00226 template <class TPManager>
00227 void
00228 QCMndo::deriveElecEnergy (TPManager& manager,
00229                           typename TPManager::TSystem::QCIterator& workingSystem,
00230                           const QCPoint3D& shift,
00231                           QCFloat * QCRestrict dRepIntegralSSVector) {
00232 
00233 
00234 #ifdef VERBOSE_DERIV2
00235   cout << " *   Derive elec energy  * " << endl;
00236 #endif
00237 
00238   // On recupere avec des restricts pour le calcul les
00239   // structures qui contiennent les donnees a manipuler
00240   QCSystem&        QCRestrict currentSystem   = workingSystem->getSystem();
00241   QCModelMatrices& QCRestrict currentMatrices = workingSystem->getMatrices();
00242   
00243   // Les parametres
00244   const typename TPManager::TParam *QCRestrict parameters = 
00245     manager.getParameters();
00246   
00247   // Le modele
00248   //   typename TPManager::TModel& QCRestrict model =
00249   //     manager.getModel();
00250 
00251   // General data
00252   const QCGeneralData& data = manager.getGeneralData();
00253 
00254   // Les atomes
00255   QCPoint3D  coordsA;
00256   QCPoint3D  coordsB;
00257 
00258   // leurs parametres
00259   const typename TPManager::TParam *QCRestrict parameterA;
00260   const typename TPManager::TParam *QCRestrict parameterB;
00261 
00262   // La matrice densite
00263   QCSymMatrix& QCRestrict densityMatrix = 
00264     currentMatrices.getDensityP().getMatrix();
00265 
00266   // La matrice des derivees des integrales d overlap.
00267   QCSymMatrix& QCRestrict dOverlapMatrix = 
00268     currentMatrices.getDOverlapSOnDr()->getMatrix();
00269 
00270   // La matrice des derivees des integrales de repulsion
00271   QCRepInterElec& QCRestrict repInterElecIntegrals =
00272     currentMatrices.getRepInterElecIntegrals();
00273 
00274   QCRepInterElec * QCRestrict dRepInterElecIntgOnDCoord =
00275     currentMatrices.getDRepInterElecIntegrals();
00276 
00277   QCMatrix& QCRestrict dRepIntegralsMatrix =
00278     currentMatrices.getDRepInterElecIntegrals()->getMatrix();
00279 
00280 
00281   //les sous-matrices
00282   QCMatrix& QCRestrict          dOverlapSAB = spdWorkingAB;
00283   QCSymMatrix& QCRestrict       densityPAA  = spWorkingAA;
00284   QCMatrix& QCRestrict          densityPAB  = spWorkingAB;
00285   QCSymMatrix& QCRestrict       densityPBB  = spWorkingBB;
00286 
00287 
00288 
00289   // distance entre a et b.
00290   QCFloat rAB;
00291 
00292   //Les derivees par rapport a r, x, y et z.
00293   QCFloat       deOnDCoord;
00294 
00295   //Les premiers index d orbitales de A et B dans la matrice complete du
00296   //systeme et l index du premier index du successeur de A tjs dans la 
00297   //matrice complete.
00298   int   firstAOofA, firstAOofB;
00299   const int SSRepIndex = 0;
00300   int i, j, mu, nu, lda, sma, muNu, ldaSma;
00301   int dRepIntegralSSIndex;
00302   register QCFloat interVal, interVal2;
00303   register QCFloat interValDiag, interValHorsDiag;
00304   register QCFloat interValDiag2, interValHorsDiag2;
00305   register QCFloat interValHorsDiag3, interValDiag3;
00306 
00307 
00308   // La double boucle sur les atomes.
00309   // je prends i < j.
00310   dRepIntegralSSIndex = 0;
00311   for (i = 1; i < currentSystem.getNbAtoms(); ++i) {
00312 
00313     currentSystem.getPointAt(i, coordsA);
00314     parameterA = &parameters[currentSystem.getParamIndexAt(i)];
00315     firstAOofA = currentSystem.getFirstAOAt(i);
00316 
00317     // Extraction de sous matrice
00318     densityMatrix.extractSubTriangle(densityPAA,
00319                                      firstAOofA, parameterA->getNbAO());
00320 
00321     // Le cas ou A est different de B.
00322     // On veut assigner les F_mu_lambda.
00323     for (j = 0; j < i; ++j) { 
00324       
00325       currentSystem.getPointAt(j, coordsB);
00326       parameterB = &parameters[currentSystem.getParamIndexAt(j)];
00327       firstAOofB = currentSystem.getFirstAOAt(j);
00328       
00329       rAB     = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00330                      QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00331                      QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00332       
00333       // Extraction de sous matrices
00334       densityMatrix.extractSubMatrix(densityPAB,
00335                                      firstAOofA, parameterA->getNbAO(),
00336                                      firstAOofB, parameterB->getNbAO() );
00337       densityMatrix.extractSubTriangle(densityPBB,
00338                                        firstAOofB, parameterB->getNbAO() );
00339 
00340 
00341       //------------------------------------------------------------//
00342       //  Partie calcul de Hamilton I
00343       //------------------------------------------------------------//
00344 
00345       // On extrait la derivee de la sous-matrice d overlap dont on a besoin.
00346       dOverlapMatrix.extractSubMatrix(dOverlapSAB,
00347                                       firstAOofA, parameterA->getNbAO(),
00348                                       firstAOofB, parameterB->getNbAO() );
00349  
00350       // Calcul de la derivee de l energie par rapport a rAB
00351       deOnDCoord = QC_ZERO;
00352 
00353       //derivee de H_AB
00354       interVal = QC_ZERO;
00355       for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00356         for (lda = 0;  lda < parameterB->getNbAO(); ++lda) {
00357           interVal += densityPAB[mu][lda] *
00358             (parameterA->getBeta(mu) + parameterB->getBeta(lda) ) * dOverlapSAB[mu][lda];
00359         }
00360       }
00361 
00362       // Attention, on ne divise pas par 2 car il existe 2 termes symetriques
00363       // qui apportent leur contribution.
00364       deOnDCoord += interVal;
00365 
00366       dRepInterElecIntgOnDCoord->computeDerivatives(parameterA,
00367                                                     parameterB,
00368                                                     repInterElecIntegrals,
00369                                                     data.getDerivationType(),
00370                                                     coordsA,
00371                                                     coordsB,
00372                                                     shift);
00373 
00374 
00375       // On utilisera ce vecteur dans deriveCoreEnergy.
00376       // Alors attention, ici on a d<ss|ss>/dx ou d<ss|ss>/dy ou d<ss|ss>/dz.
00377       // et on cherche a avoir d<ss|ss>/dr = d<ss|ss>/dx * deltaX / r + ...
00378       // soit le produit scalaire df.r/||r||, soit la projection du vecteur 
00379       // df(df/dx, df/dy, df/dz) dans la direction du vecteur r.
00380       if (shift[COORDX] != QC_ZERO) {
00381         dRepIntegralSSVector[dRepIntegralSSIndex++] += 
00382           dRepIntegralsMatrix[SSRepIndex][SSRepIndex] * 
00383           ( ( coordsB[COORDX] - coordsA[COORDX] ) / rAB );
00384       
00385       } else if (shift[COORDY] != QC_ZERO) {
00386         dRepIntegralSSVector[dRepIntegralSSIndex++] +=
00387           dRepIntegralsMatrix[SSRepIndex][SSRepIndex] * 
00388           ( ( coordsB[COORDY] - coordsA[COORDY] ) / rAB );
00389       
00390       } else if (shift[COORDZ] != QC_ZERO) {
00391         dRepIntegralSSVector[dRepIntegralSSIndex++] +=
00392           dRepIntegralsMatrix[SSRepIndex][SSRepIndex] *
00393           ( ( coordsB[COORDZ] - coordsA[COORDZ] ) / rAB );
00394       }
00395 
00396 
00397 
00398       //------------------------------------------------------------//
00399       //  Partie calcul de Hamilton II
00400       //------------------------------------------------------------//
00401       
00402       // derivee de H_BB : termes hors diag
00403       interValHorsDiag = QC_ZERO;
00404       interValDiag = QC_ZERO;
00405       for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00406         muNu = mu * parameterA->getNbAO();
00407         for (nu = 0;  nu < mu; ++nu, ++muNu) {
00408           interValHorsDiag += densityPAA[mu][nu] * dRepIntegralsMatrix[muNu][SSRepIndex];
00409         }
00410 
00411 
00412         // traitement du terme diag.
00413         // Attention, ici, je profite du fait que nu == mu en sortie
00414         // de la bcle precedente, idem pour nuNu qui est a une valeur adequate.
00415         interValDiag += densityPAA[mu][mu] * dRepIntegralsMatrix[muNu][SSRepIndex];
00416       }
00417 
00418       // mult par 2 pour les termes hors diag a cause de leur symetrique.
00419       deOnDCoord -= QC_TWO * parameterB->getCoreCharge() * interValHorsDiag;
00420       deOnDCoord -= parameterB->getCoreCharge() * interValDiag;
00421 
00422       // derivee de H_AA : termes hors diag
00423       interValHorsDiag = QC_ZERO;
00424       interValDiag = QC_ZERO;
00425       for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00426         ldaSma = lda * parameterB->getNbAO();
00427         for (sma = 0; sma < lda; ++sma, ++ldaSma) {
00428           interValHorsDiag += densityPBB[lda][sma] * dRepIntegralsMatrix[SSRepIndex][ldaSma];
00429         }
00430 
00431         // traitement du terme diag.
00432         // Attention, ici, je profite du fait que lda == sma en sortie
00433         // de la bcle precedente, idem pour ldaSma qui est a une valeur adequate.
00434         interValDiag += densityPBB[lda][lda] * dRepIntegralsMatrix[SSRepIndex][ldaSma];
00435       }
00436 
00437       // mult par 2 pour les termes hors diag a cause de leurs symetriques.
00438       deOnDCoord -= QC_TWO * parameterA->getCoreCharge() * interValHorsDiag;
00439       deOnDCoord -= parameterA->getCoreCharge() * interValDiag;
00440 
00441 
00442 
00443       //------------------------------------------------------------//
00444       //  Partie calcul de Fock
00445       //------------------------------------------------------------//
00446 
00447       interVal = QC_ZERO;
00448 
00449       // Je m occupe d abord des termes batis a partir de PAB.
00450       for (mu = 0, muNu = 0; mu < parameterA->getNbAO(); ++mu) {
00451         for (lda = 0, ldaSma = 0; lda < parameterB->getNbAO(); ++lda) {
00452           interVal2 = QC_ZERO;
00453           muNu = mu * parameterA->getNbAO();
00454           for (nu = 0;  nu < parameterA->getNbAO(); ++nu, ++muNu) {
00455             ldaSma = lda * parameterB->getNbAO();
00456             for (sma = 0; sma < parameterB->getNbAO(); ++sma, ++ldaSma) {
00457               interVal2 += densityPAB[nu][sma] * dRepIntegralsMatrix[muNu][ldaSma];
00458             }
00459           }
00460           interVal += densityPAB[mu][lda] * interVal2;
00461         }
00462       }
00463       deOnDCoord += -QC_HALF * interVal;
00464       
00465       // puis des termes batis a partir de PAA et PBB.
00466       // D abord les termes Hors Diag sur les 2 sous-triangles.
00467       interValHorsDiag = QC_ZERO;
00468       interValDiag = QC_ZERO;
00469       interValHorsDiag3 = QC_ZERO;
00470       interValDiag3 = QC_ZERO;
00471       for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00472         muNu = mu * parameterA->getNbAO();
00473         for (nu = 0;  nu < mu; ++nu, ++muNu) {
00474           interValHorsDiag2 = QC_ZERO;
00475           interValDiag2 = QC_ZERO;
00476           for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00477             ldaSma = lda * parameterB->getNbAO();
00478             for (sma = 0; sma < lda; ++sma, ++ldaSma) {
00479               interValHorsDiag2 += densityPBB[lda][sma] * dRepIntegralsMatrix[muNu][ldaSma];
00480             }
00481             interValDiag2 += densityPBB[lda][lda] * dRepIntegralsMatrix[muNu][ldaSma];
00482           }
00483           interValHorsDiag += densityPAA[mu][nu] * interValHorsDiag2;
00484           interValDiag  += densityPAA[mu][nu] * interValDiag2;
00485         }
00486 
00487         // Puis, les termes diag de PAA.
00488         interValHorsDiag2 = QC_ZERO;
00489         interValDiag2 = QC_ZERO;
00490         for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00491           ldaSma = lda * parameterB->getNbAO();
00492           for (sma = 0; sma < lda; ++sma, ++ldaSma) {
00493 
00494             // Avec les termes hors diag de PBB
00495             interValHorsDiag2 += densityPBB[lda][sma] * dRepIntegralsMatrix[muNu][ldaSma];
00496           }
00497 
00498           // Puis avec les termes diag de PBB.
00499           interValDiag2 += densityPBB[lda][lda] * dRepIntegralsMatrix[muNu][ldaSma];
00500         }
00501 
00502         // Traitement des terme diag de PAA avec les termes hors diag de PBB
00503         interValHorsDiag3 += densityPAA[mu][mu] * interValHorsDiag2;
00504 
00505         // Traitement du terme diag de PAA avec le terme diag de PBB
00506         interValDiag3 += densityPAA[mu][mu] * interValDiag2;
00507       }
00508 
00509       // On multiplie par 4 car on compte les termes doublement symetriques 4 fois
00510       deOnDCoord += QC_FOUR * interValHorsDiag;
00511       
00512       // Traitement du terme hors diag de PAA avec le terme diag de PBB
00513       deOnDCoord += QC_TWO * interValDiag;
00514 
00515       // Traitement des terme diag de PAA avec les termes hors diag de PBB
00516       deOnDCoord += QC_TWO * interValHorsDiag3;
00517 
00518       // Traitement du terme diag de PAA avec le terme diag de PBB
00519       deOnDCoord += interValDiag3;
00520 
00521 
00522       // Dispatch de cette derivee par rapport aux coordonnees
00523       dispatchDerivative (currentSystem,
00524                           (shift[COORDX] != QC_ZERO)? deOnDCoord : QC_ZERO,
00525                           (shift[COORDY] != QC_ZERO)? deOnDCoord : QC_ZERO,
00526                           (shift[COORDZ] != QC_ZERO)? deOnDCoord : QC_ZERO,
00527                           i,
00528                           j);
00529       
00530     }
00531   }
00532 
00533 
00534 }
00535 
00536 
00537 
00538 
00542 template <class TPManager>
00543 void
00544 QCMndo::deriveCoreEnergy (TPManager& manager,
00545                           typename TPManager::TSystem::QCIterator& workingSystem,
00546                           QCFloat * QCRestrict dRepIntegralSSVector) {
00547 
00548   
00549 #ifdef VERBOSE_DERIV2
00550   cout << " *   Derive core energy  * " << endl;
00551 #endif
00552 
00553   // Le modele
00554   typename TPManager::TModel& QCRestrict model =
00555     manager.getModel();
00556 
00557   // On recupere avec des restricts pour le calcul les
00558   // structures qui contiennent les donnees a manipuler
00559   QCSystem&        QCRestrict currentSystem   = workingSystem->getSystem();
00560   QCModelMatrices& QCRestrict currentMatrices = workingSystem->getMatrices();
00561 
00562   // La matrice des derivees des integrales de repulsion
00563   QCRepInterElec& QCRestrict repInterElecIntegrals =
00564     currentMatrices.getRepInterElecIntegrals();
00565   QCMatrix& QCRestrict repIntegralsMatrix =
00566     repInterElecIntegrals.getMatrix();
00567   
00568   // Les parametres
00569   const typename TPManager::TParam *QCRestrict parameters = 
00570     manager.getParameters();
00571   
00572   // General data
00573   const QCGeneralData& data = manager.getGeneralData();
00574 
00575   // Les atomes
00576   QCPoint3D  coordsA;
00577   QCPoint3D  coordsB;
00578 
00579   // leurs parametres
00580   const typename TPManager::TParam *QCRestrict parameterA;
00581   const typename TPManager::TParam *QCRestrict parameterB;
00582 
00583 
00584   // L'index de ss dans le tableau des integrales de rep inter elec.
00585   const int     SSRepIndex = 0;
00586 
00587   // Les types d atomes
00588   int atomicNbA;
00589   int atomicNbB;
00590 
00591   
00592   //distance entre a et b et son inverse.
00593   QCFloat       rAB, invrAB;
00594   QCFloat       deOnDr;
00595   QCFloat       ZAZB;
00596   QCFloat       alfAR, alfBR;
00597   QCFloat       expMalfAR, expMalfBR;
00598   QCFloat       deltaXOnR, deltaYOnR, deltaZOnR;
00599   QCFloat       dRepIntegralSS;
00600   QCFloat       inv2DeltaQ = QC_ONE / (QC_TWO * data.getDeltaQ() );
00601   int dRepIntegralSSIndex;
00602   int i, j;
00603 
00604 
00605   // La double boucle sur les atomes.
00606   // je prends i < j.
00607   deOnDr = QC_ZERO;
00608   dRepIntegralSSIndex = 0;
00609   for (i = 1; i < currentSystem.getNbAtoms(); ++i) {
00610 
00611     currentSystem.getPointAt(i, coordsA);
00612     parameterA = &parameters[currentSystem.getParamIndexAt(i)];
00613     atomicNbA  = parameterA->getAtomicNbZ();
00614 
00615     for (j = 0; j < i; ++j) {
00616 
00617       currentSystem.getPointAt(j, coordsB);
00618       parameterB = &parameters[currentSystem.getParamIndexAt(j)];
00619       atomicNbB  = parameterB->getAtomicNbZ();
00620       
00621       rAB     = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00622                      QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00623                      QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00624 
00625       invrAB = QC_ONE / rAB;
00626 
00627       deltaXOnR = (coordsB[COORDX]-coordsA[COORDX]) * invrAB;
00628       deltaYOnR = (coordsB[COORDY]-coordsA[COORDY]) * invrAB;
00629       deltaZOnR = (coordsB[COORDZ]-coordsA[COORDZ]) * invrAB;
00630 
00631       dRepIntegralSS = dRepIntegralSSVector[dRepIntegralSSIndex++] * inv2DeltaQ;
00632 
00633       //Calcul des elements de la matrice de rep inter elec.
00634       repInterElecIntegrals.computeElems(parameterA,
00635                                          parameterB,
00636                                          coordsA,
00637                                          coordsB,
00638                                          rAB,
00639                                          //
00640                                          //ssACalc recoit true,
00641                                          true,
00642                                          //
00643                                          //ssBCalc recoit true,
00644                                          //car ici, on ne fait les calculs que pour 
00645                                          //l element <ss|ss> de la matrice de rep
00646                                          //inter-elec.
00647                                          true);
00648 
00649 
00650       //   
00651       ZAZB      = parameterA->getCoreCharge() * parameterB->getCoreCharge();
00652       alfAR     = parameterA->getAlpha() * rAB;
00653       expMalfAR = QCExp(-alfAR);
00654       alfBR     = parameterB->getAlpha() * rAB;
00655       expMalfBR = QCExp(-alfBR);
00656 
00657 
00658       // Cas ou on n a pas une liaison N-H ou O-H.
00659       // Attention, ici, pour le test, on ne peut pas se servir de 
00660       // sometimesR a cause du cas ambigu ou rAB == QCOne.
00661       if (isOHorNH(atomicNbA, atomicNbB) ) {
00662         deOnDr = ZAZB *
00663           (dRepIntegralSS * (QC_ONE + rAB * expMalfAR + expMalfBR) +
00664            repIntegralsMatrix[SSRepIndex][SSRepIndex] * ( (QC_ONE - alfAR) * expMalfAR +
00665                                                           (-parameterB->getAlpha() ) * expMalfBR) );
00666       }
00667       else if (isOHorNH(atomicNbB, atomicNbA) ) {
00668         deOnDr = ZAZB *
00669           (dRepIntegralSS * (QC_ONE + expMalfAR + rAB * expMalfBR) +
00670            repIntegralsMatrix[SSRepIndex][SSRepIndex] * ( (-parameterA->getAlpha() ) * expMalfAR +
00671                                                           (QC_ONE - alfBR) * expMalfBR) );
00672       }
00673       else {//on n a pas une liaison O-H ou N-H.
00674         deOnDr = ZAZB *
00675           (dRepIntegralSS * (QC_ONE + expMalfAR + expMalfBR) +
00676            repIntegralsMatrix[SSRepIndex][SSRepIndex] * ( (-parameterA->getAlpha() ) * expMalfAR +
00677                                                           (-parameterB->getAlpha() ) * expMalfBR) );
00678       }
00679 
00680       deOnDr += 
00681         (ZAZB * invrAB) *
00682         (-model.exponentialExpansion(parameterA, parameterB, rAB) * invrAB +
00683          model.dExponentialExpansion(parameterA, parameterB, rAB) );
00684 
00685       // Dispatch de cette derivee par rapport a en derivees par rapport aux coordonnees.
00686       // Attention, on forwarde de/dxj = de/dr * (xj-xi)/r ainsi que de/dyj et de/dzj.
00687       dispatchDerivative(currentSystem,
00688                          deOnDr * deltaXOnR,
00689                          deOnDr * deltaYOnR,
00690                          deOnDr * deltaZOnR,
00691                          i, j);
00692 
00693     }
00694   }
00695 
00696 
00697 
00698 }
00699 
00700 
00701 
00702 
00703 
00707 void
00708 QCMndo::dispatchDerivative (QCSystem&  system,
00709                             QCFloat    deOnDxj,
00710                             QCFloat    deOnDyj,
00711                             QCFloat    deOnDzj,
00712                             int        i, 
00713                             int        j) {
00714 
00715   if (deOnDxj != QC_ZERO) {
00716     system.incrXDisp(i, -deOnDxj);
00717     system.incrXDisp(j, deOnDxj);
00718   }
00719   if (deOnDyj != QC_ZERO) {
00720     system.incrYDisp(i, -deOnDyj);
00721     system.incrYDisp(j, deOnDyj);
00722   }
00723   if (deOnDzj != QC_ZERO) { 
00724     system.incrZDisp(i, -deOnDzj);
00725     system.incrZDisp(j, deOnDzj);
00726   }
00727 }
00728 
00729 
00730 
00731 
00732 
00736 template <class TPParam>
00737 QCMndoT<TPParam>::QCMndoT (void) :
00738   QCMndo( (TPParam::NB_S_ORB + TPParam::NB_P_ORB),
00739           (TPParam::NB_S_ORB + TPParam::NB_P_ORB + TPParam::NB_D_ORB) )
00740 {}
00741 
00742 
00743 
00744 
00748 template <class TPParam>
00749 QCMndoT<TPParam>::~QCMndoT () 
00750 {}
00751 
00752 
00753 
00754 
00758 template class QCMndoT<QCMndoParam>;
00759 
00760 QCMANAGER_ITER_METH_EXPL_INST(QCFloat QCMndo::computeCoreEnergy);
00761 
00762 QCITERATOR_METH_EXPL_INST_PARAM(QCFloat QCMndo::computeElecEnergy,
00763                                 QCDensityGtr<QCSymMatrix>& QCRestrict);
00764 
00765 
00766 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCMndo::deriveElecEnergy,
00767                                     TWO_PARAMS(const QCPoint3D&,
00768                                                QCFloat * QCRestrict));
00769 
00770 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCMndo::deriveCoreEnergy,
00771                                     QCFloat * QCRestrict);

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