QCRepInterElec.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 //#include <stdio.h>
00019 //
00020 #include "QCRepInterElec.hpp"
00021 #include "QCIntgReader.hpp"
00022 #include "QCIntgWriter.hpp"
00023 #include "QCFortran.hpp"
00024 #include "QCMndoParam.hpp"
00025 #include "QCAm1Param.hpp"
00026 #include "QCPm3Param.hpp"
00027 #include "QCSystem.hpp"
00028 
00029 
00030 
00034 enum QCRepInterElecInd10 {
00035   s_s_10,       s_px_10,        s_py_10,        s_pz_10,
00036   px_px_10,     px_py_10,       px_pz_10,
00037   py_py_10,     py_pz_10,
00038   pz_pz_10
00039 };
00040 enum QCRepInterElecInd16 {
00041   s_s_16,       s_px_16,        s_py_16,        s_pz_16,
00042   px_s_16,      px_px_16,       px_py_16,       px_pz_16,
00043   py_s_16,      py_px_16,       py_py_16,       py_pz_16,  
00044   pz_s_16,      pz_px_16,       pz_py_16,       pz_pz_16
00045 };
00046 
00047 
00048 
00049 
00053 QCRepInterElec::QCRepInterElec (const int& dim) : 
00054   QCMatElemGtr<QCMatrix>(dim),
00055   intgWriter(NULL),
00056   intgReader(NULL),
00057   intgTmpFile("") 
00058 {}
00059 
00060 
00061 
00065 QCRepInterElec::~QCRepInterElec (void) {
00066   if (intgWriter) { delete intgWriter; intgWriter = NULL; }
00067   if (intgReader) { delete intgReader; intgReader = NULL; }
00068 
00069   if (intgTmpFile != "") {
00070     //
00071     // Remove tempory file (STDlib) 
00072     remove(intgTmpFile.c_str());
00073   }
00074 }
00075 
00076 
00077 
00081 QCIntgWriter * 
00082 QCRepInterElec::allocIntgWriter (void) {
00083   if (!intgWriter) {
00084     intgWriter = new QCIntgWriter;
00085   }
00086   return intgWriter; 
00087 }
00088 
00089 
00090 
00094 QCIntgReader * 
00095 QCRepInterElec::allocIntgReader (void) {
00096   if (!intgReader) {
00097     intgReader = new QCIntgReader;
00098   }
00099   return intgReader;
00100 }
00101 
00102 
00103 
00108 void
00109 QCRepInterElec::expandOverDim16 (const int& nbAOofA, const int& nbAOofB) {
00110   //
00111   //l index de ss dans le tableau des integrales de rep inter elec.
00112   const int SSRepIndex    = 0;
00113   const int compressedDim = 10;
00114   
00115   if (nbAOofA == 1) {
00124     matrix[SSRepIndex][pz_pz_16] = matrix[SSRepIndex][pz_pz_10];
00125     matrix[SSRepIndex][pz_py_16] = matrix[SSRepIndex][py_pz_10];
00126     matrix[SSRepIndex][pz_px_16] = matrix[SSRepIndex][px_pz_10];
00127     matrix[SSRepIndex][pz_s_16]  = matrix[SSRepIndex][s_pz_10];
00133     matrix[SSRepIndex][py_py_16] = matrix[SSRepIndex][py_py_10];
00134     matrix[SSRepIndex][px_py_16] = matrix[SSRepIndex][px_py_10];
00135     matrix[SSRepIndex][py_px_16] = matrix[SSRepIndex][px_py_10];
00136     matrix[SSRepIndex][px_px_16] = matrix[SSRepIndex][px_px_10];
00140     matrix[SSRepIndex][px_s_16]  = matrix[SSRepIndex][s_px_16];
00141     matrix[SSRepIndex][py_s_16]  = matrix[SSRepIndex][s_py_16];
00142     matrix[SSRepIndex][px_pz_16] = matrix[SSRepIndex][pz_px_16];
00143     matrix[SSRepIndex][py_pz_16] = matrix[SSRepIndex][pz_py_16];
00144   
00145   } else if (nbAOofB == 1) {
00149     matrix[pz_pz_16][SSRepIndex] = matrix[pz_pz_10][SSRepIndex];
00150     matrix[pz_py_16][SSRepIndex] = matrix[py_pz_10][SSRepIndex];
00151     matrix[pz_px_16][SSRepIndex] = matrix[px_pz_10][SSRepIndex];
00152     matrix[pz_s_16][SSRepIndex]  = matrix[s_pz_10][SSRepIndex];
00158     matrix[py_py_16][SSRepIndex] = matrix[py_py_10][SSRepIndex];
00159     matrix[px_py_16][SSRepIndex] = matrix[px_py_10][SSRepIndex];
00160     matrix[py_px_16][SSRepIndex] = matrix[px_py_10][SSRepIndex];
00161     matrix[px_px_16][SSRepIndex] = matrix[px_px_10][SSRepIndex];
00165     matrix[px_s_16][SSRepIndex]  = matrix[s_px_16][SSRepIndex];
00166     matrix[py_s_16][SSRepIndex]  = matrix[s_py_16][SSRepIndex];
00167     matrix[px_pz_16][SSRepIndex] = matrix[pz_px_16][SSRepIndex];
00168     matrix[py_pz_16][SSRepIndex] = matrix[pz_py_16][SSRepIndex]; 
00169   
00170   } else { // nbAOofA != 1 && nbAOofB != 1
00174     for (int i = 0; i < compressedDim; ++i) {
00181       matrix[i][pz_pz_16] = matrix[i][pz_pz_10];
00182       matrix[i][pz_py_16] = matrix[i][py_pz_10];
00183       matrix[i][pz_px_16] = matrix[i][px_pz_10];
00184       matrix[i][pz_s_16]  = matrix[i][s_pz_10];
00190       matrix[i][py_py_16] = matrix[i][py_py_10];
00191       matrix[i][px_py_16] = matrix[i][px_py_10];
00192       matrix[i][py_px_16] = matrix[i][px_py_10];
00193       matrix[i][px_px_16] = matrix[i][px_px_10];
00197       matrix[i][px_s_16]  = matrix[i][s_px_16];
00198       matrix[i][py_s_16]  = matrix[i][s_py_16];
00199       matrix[i][px_pz_16] = matrix[i][pz_px_16];
00200       matrix[i][py_pz_16] = matrix[i][pz_py_16];
00201     }
00211     memcpy(matrix[pz_pz_16], matrix[pz_pz_10], matrix.getDim()*sizeof(QCFloat) );
00212     memcpy(matrix[pz_py_16], matrix[py_pz_10], matrix.getDim()*sizeof(QCFloat) );
00213     memcpy(matrix[pz_px_16], matrix[px_pz_10], matrix.getDim()*sizeof(QCFloat) );
00214     memcpy(matrix[pz_s_16],  matrix[s_pz_10],  matrix.getDim()*sizeof(QCFloat) );
00220     memcpy(matrix[py_py_16], matrix[py_py_10], matrix.getDim()*sizeof(QCFloat) );  
00221     memcpy(matrix[px_py_16], matrix[px_py_10], matrix.getDim()*sizeof(QCFloat) );
00222     memcpy(matrix[py_px_16], matrix[px_py_10], matrix.getDim()*sizeof(QCFloat) );
00223     memcpy(matrix[px_px_16], matrix[px_px_10], matrix.getDim()*sizeof(QCFloat) );
00227     memcpy(matrix[px_s_16],  matrix[s_px_16],  matrix.getDim()*sizeof(QCFloat) );  
00228     memcpy(matrix[py_s_16],  matrix[s_py_16],  matrix.getDim()*sizeof(QCFloat) );
00229     memcpy(matrix[px_pz_16], matrix[pz_px_16], matrix.getDim()*sizeof(QCFloat) );
00230     memcpy(matrix[py_pz_16], matrix[pz_py_16], matrix.getDim()*sizeof(QCFloat) );  
00231   }
00232 }
00233 
00234 
00235 
00236 
00241 template <class TPParam>
00242 void 
00243 QCRepInterElec::computeElems (const TPParam*   QCRestrict parameterA,
00244                               const TPParam*   QCRestrict parameterB,
00245                               const QCPoint3D& QCRestrict coordsA,
00246                               const QCPoint3D& QCRestrict coordsB,
00247                               QCFloat rAB,
00248                               bool  ssACalc,
00249                               bool  ssBCalc,
00250                               bool  compressedFormat,
00251                               const QCIntgAcquisitionMethod& intgAcquisitionMethod,
00252                               bool  isItfirstCall)
00253 {
00254   
00255   int nbComputedAOofA;
00256   int nbComputedAOofB;
00257   int compressedFormatInt;
00262   if (compressedFormat) {
00263     compressedFormatInt = 1;
00264   
00265   } else {
00266     compressedFormatInt = 0;
00267   }
00268   
00269   if (ssBCalc) {
00270     nbComputedAOofB = TPParam::getNbSOrb();
00271   
00272   } else {
00273     nbComputedAOofB = parameterB->getNbAO();
00274   }
00275   
00284   if (ssACalc) {
00285     nbComputedAOofA = TPParam::getNbSOrb();
00286   
00287   } else {
00288     nbComputedAOofA = parameterA->getNbAO();
00289   }
00290 
00297   if ( (intgAcquisitionMethod == QC_DIRECT_COMPUTATION) ||
00298        (intgAcquisitionMethod == QC_INDIRECT_STORAGE && isItfirstCall) ) {
00299     
00300     QCFortran(repul)
00301       (compressedFormatInt,
00302        // atom B
00303        parameterB->getAtomicNbZ(),
00304        //Attention, pour B, on ne fait le calcul
00305        //que sur la rangee ss.
00306        nbComputedAOofB,
00307        parameterB->getRo0(),
00308        parameterB->getRo1(),
00309        parameterB->getRo2(),
00310        parameterB->getD1(),
00311        parameterB->getD2(),
00312        coordsB[COORDX],
00313        coordsB[COORDY],
00314        coordsB[COORDZ],
00315        // atom A
00316        parameterA->getAtomicNbZ(),
00317        nbComputedAOofA,
00318        parameterA->getRo0(),
00319        parameterA->getRo1(),
00320        parameterA->getRo2(),
00321        parameterA->getD1(),
00322        parameterA->getD2(),
00323        coordsA[COORDX],
00324        coordsA[COORDY],
00325        coordsA[COORDZ],
00326        rAB,
00327        matrix.getElems() );
00328     
00329     if (intgAcquisitionMethod == QC_INDIRECT_STORAGE) {
00330 
00331       // Attention, on compresse pour avoir moins a stocker.
00332       if (nbComputedAOofA > nbComputedAOofB) {
00333         matrix.pack(QCPow<2>(nbComputedAOofA) );
00334       }
00335 
00336       //On ecrit dans le fichier de stockage.
00337       intgWriter->storeIntegrals(matrix.getElems(), 
00338                                  sizeof(QCFloat),
00339                                  QCPow<2>(nbComputedAOofA) * QCPow<2>(nbComputedAOofB) );
00340       
00341       if (nbComputedAOofA > nbComputedAOofB) {
00342         matrix.unpack(QCPow<2>(nbComputedAOofA) );
00343       }
00344     }
00345     
00346   } else { 
00347     //(intgAcquisitionMethod == QC_DIRECT_COMPUTATION) ||
00348     //(intgAcquisitionMethod == QC_INDIRECT_STORAGE && isItfirstCall)
00349 
00350     // Donc on lit.
00351     intgReader->loadIntegrals(matrix.getElems(), 
00352                               sizeof(QCFloat),
00353                               QCPow<2>(nbComputedAOofA) * QCPow<2>(nbComputedAOofB) );
00354 
00355     // Ne pas oublier de decompresser.
00356     if (nbComputedAOofA > nbComputedAOofB) {
00357       matrix.unpack(QCPow<2>(nbComputedAOofA) );
00358     }
00359   }  
00360 
00361   
00362 }
00363 
00364 
00365 
00366 
00367 
00371 template <class TPParam>
00372 void  
00373 QCRepInterElec::computeDerivatives (const TPParam*    QCRestrict parameterA,
00374                                     const TPParam*    QCRestrict parameterB, 
00375                                     QCRepInterElec&   QCRestrict workingIntg,
00376                                     QCDerivationType  derivationType,
00377                                     const QCPoint3D&  QCRestrict coordsA,
00378                                     const QCPoint3D&  QCRestrict coordsB,
00379                                     const QCPoint3D&  shift,
00380                                     bool              ssACalc,
00381                                     bool              ssBCalc,
00382                                     const QCFloat&    QCRestrict inv2DeltaQ
00383 #ifdef TRACE_COMPUTE_DERIVATIVES
00384                                     ,
00385                                     QCChrono *        intgChr1Ptr = NULL,
00386                                     QCChrono *        intgChr2Ptr = NULL, 
00387                                     QCChrono *        intgChr3Ptr = NULL
00388 #endif // TRACE_COMPUTE_DERIVATIVES
00389                                     ) {
00390   
00391   int       nbAOofA, nbAOofB;
00392   QCPoint3D coordsBdelta;
00393   QCFloat   rAB;
00394   const int SSRepIndex = 0;
00401   bool compressedFormat = true;
00402   
00403   if (derivationType == QC_NUMERICAL_DERIVATION) {
00404     
00405     nbAOofA = parameterA->getNbAO();
00406     nbAOofB = parameterB->getNbAO();
00407 
00408     coordsBdelta  = coordsB;
00409     coordsBdelta -= shift;
00410 
00411     rAB     = sqrt(QCPow<2>(coordsBdelta[COORDX] - coordsA[COORDX]) +
00412                    QCPow<2>(coordsBdelta[COORDY] - coordsA[COORDY]) +
00413                    QCPow<2>(coordsBdelta[COORDZ] - coordsA[COORDZ]) );
00414     
00415 #ifdef TRACE_COMPUTE_DERIVATIVES
00416     if (intgChr1Ptr) { intgChr1Ptr->run(); }
00417 #endif // TRACE_COMPUTE_DERIVATIVES
00418 
00419     workingIntg.computeElems(// Attention, le premier entree
00420                              // est celui dont l indice varie le plus vite.
00421                              parameterA,
00422                              parameterB,
00423                              coordsA,
00424                              coordsBdelta,
00425                              rAB,
00426                              // ssACalc et ssBCalc recoivent false,
00427                              // car ici, on fait tous les calculs
00428                              // des elements de la matrice de rep inter elec.
00429                              ssACalc,
00430                              ssBCalc,
00431                              compressedFormat,
00432                              // calcul direct, pas de stockage disque 
00433                              // quoiqu il arrive.
00434                              QC_DIRECT_COMPUTATION);
00435 
00436     
00437 #ifdef TRACE_COMPUTE_DERIVATIVES
00438     if (intgChr1Ptr) { intgChr1Ptr->pause(); }
00439 #endif // TRACE_COMPUTE_DERIVATIVES
00440 
00441     coordsBdelta  = coordsB;
00442     coordsBdelta += shift;
00443     
00444     rAB     = sqrt(QCPow<2>(coordsBdelta[COORDX] - coordsA[COORDX]) +
00445                    QCPow<2>(coordsBdelta[COORDY] - coordsA[COORDY]) +
00446                    QCPow<2>(coordsBdelta[COORDZ] - coordsA[COORDZ]) );
00447 
00448 #ifdef TRACE_COMPUTE_DERIVATIVES
00449     if (intgChr2Ptr) { intgChr2Ptr->run(); }
00450 #endif // TRACE_COMPUTE_DERIVATIVES
00451 
00452     computeElems(// Attention, le premier entree
00453                  // est celui dont l indice varie le plus vite.
00454                  parameterA,
00455                  parameterB,
00456                  coordsA,
00457                  coordsBdelta,
00458                  rAB,
00459                  // ssACalc et ssBCalc recoivent false,
00460                  // car ici, on fait tous les calculs
00461                  // des elements de la matrice de rep inter elec.
00462                  ssACalc,
00463                  ssBCalc,
00464                  compressedFormat,
00465                  // calcul direct, pas de stockage disque 
00466                  // quoiqu il arrive.
00467                  QC_DIRECT_COMPUTATION);
00468     
00469 #ifdef TRACE_COMPUTE_DERIVATIVES
00470     if (intgChr2Ptr) { intgChr2Ptr->pause(); }
00471     if (intgChr3Ptr) { intgChr3Ptr->run(); }
00472 #endif // TRACE_COMPUTE_DERIVATIVES
00473 
00474     if (ssACalc && ssBCalc) {
00475       
00480       if (inv2DeltaQ != QC_ZERO) {
00481         matrix[SSRepIndex][SSRepIndex] = 
00482           (matrix[SSRepIndex][SSRepIndex] - workingIntg.getMatrix()[SSRepIndex][SSRepIndex]) *
00483           inv2DeltaQ;
00484       
00485       } else {
00486        
00487         matrix[SSRepIndex][SSRepIndex] = 
00488           (matrix[SSRepIndex][SSRepIndex] - workingIntg.getMatrix()[SSRepIndex][SSRepIndex]);
00489       }
00490     
00491     } else if (ssACalc && !ssBCalc) {
00492       
00493       //Et pour finir, on calcule les derivees.
00494       matrix.calculateNumericalDerivatives(workingIntg.getMatrix(), 
00495                                            inv2DeltaQ,
00496                                            1,
00497                                            (nbAOofB * (nbAOofB+1) ) / 2);
00498       expandOverDim16(1, nbAOofB);
00499     }
00500 
00501 #ifdef DEBUG_QC
00502     else if (!ssACalc && ssBCalc) {
00503      
00504       //sinon l appli est out.
00505       cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00506       abort();
00507     }
00508 #endif //DEBUG_QC
00509     
00510     else {
00511       
00512       //Et pour finir, on calcule les derivees.
00513       matrix.calculateNumericalDerivatives(workingIntg.getMatrix(), 
00514                                            inv2DeltaQ,
00515                                            (nbAOofA * (nbAOofA+1) ) / 2 ,
00516                                            (nbAOofB * (nbAOofB+1) ) / 2);
00517       expandOverDim16(nbAOofA, nbAOofB);
00518     }
00519 
00520 #ifdef TRACE_COMPUTE_DERIVATIVES
00521     if (intgChr3Ptr) { intgChr3Ptr->pause(); }
00522 #endif // TRACE_COMPUTE_DERIVATIVES
00523     
00524     
00525   } else { // derivationType == QC_ANALYTICAL_DERIVATION
00526     
00527   }
00528   
00529 }
00530 
00531 
00532 
00533 
00534 
00535 
00536 
00537 
00538 
00539 
00540 
00544 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeElems,
00545                                       EIGHT_PARAMS(const QCPoint3D& QCRestrict,
00546                                                    const QCPoint3D& QCRestrict,
00547                                                    QCFloat,
00548                                                    bool,
00549                                                    bool,
00550                                                    bool,
00551                                                    const QCIntgAcquisitionMethod&,
00552                                                    bool));
00553 
00554 #ifndef TRACE_COMPUTE_DERIVATIVES  
00555 
00556 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeDerivatives,
00557                                 EIGHT_PARAMS(QCRepInterElec&,
00558                                              QCDerivationType,
00559                                              const QCPoint3D&,
00560                                              const QCPoint3D&,
00561                                              const QCPoint3D&,
00562                                              bool,
00563                                              bool,
00564                                              const QCFloat& QCRestrict));
00565 
00566 #else // TRACE_COMPUTE_DERIVATIVES
00567 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeDerivatives,
00568                                       ELEVEN_PARAMS(QCRepInterElec&,
00569                                                     QCDerivationType,
00570                                                     const QCPoint3D&,
00571                                                     const QCPoint3D,
00572                                                     const QCPoint3D&,
00573                                                     bool,
00574                                                     bool,
00575                                                     const QCFloat&   QCRestrict 
00576                                                     QCChrono*,
00577                                                     QCChrono*,
00578                                                     QCChrono*));
00579 
00580 #endif // TRACE_COMPUTE_DERIVATIVES

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