QCSubDomain.hpp

Go to the documentation of this file.
00001 // This is a -*- C++ -*-  file
00002 //*****************************************************************************//
00003 //                                                                             //
00004 //   Copyright (c) 2001                                                        //
00005 //      INRIA                                                                  //
00006 //      54600 VILLERS LES NANCY                                                //
00007 //      France                                                                 //
00008 //                                                                             //
00009 //*****************************************************************************//
00010 //                                                                             //
00011 //               *** NOTICE OF PROPRIETARY INFORMATION ***                     //
00012 //                                                                             //
00013 // The information contained in this file is considered proprietary and the    //
00014 // exclusive property of  INRIA. This information may not be disclosed,        //
00015 // duplicated or used, in whole or in part, for  any purpose  whatsoever       //
00016 // without express written authorization from INRIA                            //
00017 //                                                                             //
00018 //*****************************************************************************//
00019 
00020 
00021 #ifndef QC_SUBDOMAIN_H_
00022 #define QC_SUBDOMAIN_H_
00023 
00024 #include <fstream>
00025 #include <vector>
00026 //
00027 #include "QCCommon.hpp"
00028 #include "QCGeneralData.hpp"
00029 #include "QCSystem.hpp"
00030 #include "QCModelMatrices.hpp"
00031 #include "QCEnergies.hpp"
00032 #include "QCPartSubDomain.hpp"
00033 
00034 
00038 #define   ATOM_OFFSET          2
00039 #define   DOM_ID_IDX           1
00040 #define   NB_ATOM_IDX          0
00041 
00042 
00043 
00047 #define   OVRLP_MAP_HEAD       1
00048 #define   OVRLP_MAP_DECL       1
00049 
00050 #define   NB_OVRLP_IDX         0
00054 #define   IND_MAP_HEAD         7
00055 #define   IND_MAP_CELL         3
00059 #define   SBD_ID_IDX           0  
00060 #define   NB_SHARED_ATOM_IDX   1  
00061 #define   NB_ATOM_SHELL1_IDX   2  
00062 #define   NB_ATOM_SHELL2_IDX   3  
00063 #define   NB_SHARED_AO_IDX     4  
00064 #define   NB_AO_SHELL1_IDX     5  
00065 #define   NB_AO_SHELL2_IDX     6  
00066 //
00067 #define   LOCAL_IND_IDX        0
00068 #define   REMOTE_IND_IDX       1
00069 #define   REMOTE_ZONE_IDX      2
00074 
00079 //                                                                  
00080 class QCSubDomain : public QCSystem, public QCModelMatrices,  
00081                     public QCEnergies
00082 {
00083   
00084   friend ostream&
00085   operator << (ostream& out, const QCSubDomain& domain);
00086 
00088 public:
00089   //-------------------------------------------------------------//
00090   // Public Method : constructor
00091   //-------------------------------------------------------------//
00092   //
00094   QCSubDomain (void) ;
00098   QCSubDomain (const QCSystem *root);
00099   //-------------------------------------------------------------//
00100   // Public Method : destructor
00101   //-------------------------------------------------------------//
00102   ~QCSubDomain  (void);
00106 
00107   int  getId (void)                 const 
00108      { return domainId; }
00111   int  getNbOverlpSD (void)         const 
00112      { return overlapMap[NB_OVRLP_IDX]; }
00114   int  getGlobalIdx (const int idx)       const { 
00115     int ret = -1;
00116     if (globalIdx) {
00117       ret = globalIdx[idx];
00118     } 
00119       return ret;
00120   }
00121   //
00122   int  getNbAtomCoreShell1 (void)   const { return nbAtomCoreShell1; }
00123   int  getNbAtomShell2 (void)       const { return nbAtomShell2; }
00124   int  getNbAOCoreShell1 (void)     const { return nbAOCoreShell1; }
00125   int  getNbAOShell2 (void)         const { return nbAOShell2; }
00126   int  getDensityBlockIdx (int idx) const { return densityBlockIdx[idx]; }
00127   int  getDensityBlockCount (void)  const { return densityBlockCount; }
00129   QCSubDomainZone getZone (int idx) const { return zone[idx]; }
00131   QCSubDomainZone* getZone ()        const { return zone; }
00132   const int * getGlobalIdxs (void)  const { return globalIdx; }
00133   QCFloat * getBFactor (void)             { return bFactor; }
00134   QCFloat   getBFactor (int idx)    const { return bFactor[idx]; }
00135   int  getNbFilledOM (void)         const { return nbFilledOM; }
00136   int  getNbOccupiedOM (void)       const { return nbOccupiedOM; }
00137   QCFloat * getOrbitalOccupN (void)       { return orbitalOccupN; }
00138   QCFloat getOrbitalOccupN (int idx) const { return orbitalOccupN[idx]; }  
00139   int  getOverlapMapSize (void)     const { return overlapMapSize; }
00140   int * getOverlapMap (void)              { return overlapMap; }
00141   
00142 
00146   void setRootSystem (const QCSystem *root)   { rootSystem = root; }
00147   void setZone (int idx, QCSubDomainZone val) { zone[idx] = val; }
00148   void setGlobalIdx (int idx, int val)        { globalIdx[idx] = val; }
00149   void setDensityBlockIdx (int idx, int val)  { densityBlockIdx[idx] = val; }
00150   void setId (int id)                         { domainId = id; }
00151   void setGlobalIndex (const int id)          { globalDomainId = id; }
00152 
00153 
00154   void setOrbitalOccupN (int idx, QCFloat val) { orbitalOccupN[idx] = val; }
00155   
00156   void setNbFilledOccupiedOM (int filled, int occup) { 
00157     nbFilledOM   = filled;   nbOccupiedOM = occup;
00158   }
00162   QCSystem& getSystem (void) { return QCSystem::asBase(); }
00166   const QCSystem& getSystem (void) const { return QCSystem::asBase(); }
00170   const QCSystem& getRootSystem (void) const { return *rootSystem; }
00174   QCModelMatrices& getMatrices (void) { return QCModelMatrices::asBase(); }
00175 
00176 
00177   void writeOverlapStructure(){
00178     std::cout <<"*******************************************" <<std::endl;
00179     for(int i =0 ; i <overlapMapSize ; ++i){std::cout << "  "<<overlapMap[i] ; }
00180     std::cout <<std::endl <<"*******************************************" <<std::endl;
00181   }
00182 
00186   int  getNbSharedAtoms (int mapidx) const {
00187     int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00188     return overlapMap[offset + NB_SHARED_ATOM_IDX];
00189   }
00190  
00194   int  getRemDomainId (int mapidx) const {
00195     int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00196     return overlapMap[offset + SBD_ID_IDX];
00197   }
00201   int  getIndLocalIdx (int mapidx, int cellidx) {
00202     int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00203     int cell   = offset + IND_MAP_HEAD + cellidx * IND_MAP_CELL;
00204     return overlapMap[cell + LOCAL_IND_IDX];
00205   }
00209   int  getIndRemoteIdx (int mapidx, int cellidx) {
00210     int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00211     int cell   = offset + IND_MAP_HEAD + cellidx * IND_MAP_CELL;
00212     return overlapMap[cell + REMOTE_IND_IDX];
00213   }
00214 
00218   int getMapOffset (int sbdId) const {
00219     int ret = -1;
00220     int offset;
00221     for (int i=0; i < getNbOverlpSD(); ++i) {
00222       offset = overlapMap[OVRLP_MAP_HEAD + i];
00223       if (overlapMap[offset + SBD_ID_IDX] == sbdId) {
00224         ret = offset;
00225         break;
00226       }
00227     }
00228     return ret;
00229   }
00230     
00234   int  getRemoteIdx (int sbdId, int atomIdx) const {
00235     int idx = -1;
00236     int offset = getMapOffset(sbdId);
00237     if (offset >= 0) {
00238       for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00239         int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00240         if (overlapMap[cellOffset + LOCAL_IND_IDX] == atomIdx) {
00241           idx = overlapMap[cellOffset + REMOTE_IND_IDX];
00242           break;
00243         }
00244       }
00245     }
00246     return idx;
00247   }
00251   bool containsNonShell2 (int mapidx) const {
00252     bool ret = false;
00253     int offset = overlapMap[OVRLP_MAP_HEAD + mapidx]; 
00254     for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00255       int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00256       if (getZone(overlapMap[cellOffset + LOCAL_IND_IDX]) != QC_SHELL2) {
00257         ret = true;
00258         break;
00259       }
00260     }
00261         
00262     return ret;
00263   }
00264   
00268   void sortOverlapments (void) {
00269     for (int i=0; i < getNbOverlpSD(); ++i) {
00270       int offset = overlapMap[OVRLP_MAP_HEAD + i];
00271       sort(offset);
00272     }
00273   }
00283   void allocStructures (int nbatoms, int nbao, int nbovrlp, const int * sbdIds, const int * nbshared);
00284 
00285   void allocStructures (QCPartSubDomain & domain );
00286     
00290   void addOverlpIndirectionInfo (int remoteSbdId,
00291                                  int atomIdx,
00292                                  int remoteAtomIdx,
00293                                  QCSubDomainZone remoteZone);
00297   template <class TPParam>
00298   void setSystemNumbers (const TPParam *QCRestrict parameters);
00302   void fillWeight (const QCGeneralData& data);
00303 
00307   template <class TPParam>
00308   void fillDensityArray (const TPParam * QCRestrict  parameters,
00309                          QCFloat  * QCRestrict &     densityArray,
00310                          bool                        alloc = false);
00314   template <class TPParam>
00315   void initWeightedDensity (const TPParam * QCRestrict);
00316 
00317 
00318 
00319 
00323   template <class TPParam>
00324   void computeBFactor (const TPParam * QCRestrict parameters);
00325 
00326 
00327 
00328 
00332   template <class TPParam>
00333   QCFloat traceProductCoreShell1 (const TPParam * QCRestrict parameters,
00334                                   QCSymMatrix&    QCRestrict fockMatrix,
00335                                   QCSymMatrix&    QCRestrict densityMatrix);
00336 
00337 
00338 
00342   void modifyMatricesAndEnergy (const QCFloat& lambdaOpt);
00346   void copyAtomOverlpInto (QCFloat * atomsBuff, int  * orderBuff,  int *ovrlpBuff,
00347                            QCSubDomainZone * zoneBuff,  int * blockBuff = NULL) {
00348 
00349     atomsBuff[NB_ATOM_IDX] = nbAtoms;
00350     atomsBuff[DOM_ID_IDX]  = domainId;
00351     for (int i=0; i < getNbAtoms() * DIMENSION; ++i) {
00352       atomsBuff[i + ATOM_OFFSET] = atoms[i];
00353     }
00354 
00355     for (int i=0; i < getNbAtoms() * ORDR_CELL_SIZE; ++i) {
00356       orderBuff[i] = cells[i];
00357     }
00358 
00359     for (int i=0; i < getOverlapMapSize(); ++i) {
00360       ovrlpBuff[i] = overlapMap[i];
00361     }
00362 
00363     for (int i=0; i < getNbAtoms(); ++i) {
00364       zoneBuff[i] = zone[i];
00365     }
00366 
00367     if (blockBuff) {
00368       for (int i=0; i < getNbAtoms(); ++i) {
00369         blockBuff[i] = densityBlockIdx[i];
00370       }
00371     }
00372   }
00376   void setAtomOverlpFrom (QCFloat * atomsBuff, int * orderBuff,  int * ovrlpBuff,
00377                           QCSubDomainZone   * zoneBuff,  int     * blockBuff = NULL) {
00378     
00379     if (atomsBuff) {
00380       nbAtoms    = static_cast<int>(atomsBuff[NB_ATOM_IDX]);
00381       domainId   = static_cast<int>(atomsBuff[DOM_ID_IDX]);
00382       atoms      = atomsBuff + ATOM_OFFSET;
00383     }
00384     cells      = orderBuff;
00385     overlapMap = ovrlpBuff;
00386     zone       = zoneBuff;
00387     if (blockBuff) {
00388       densityBlockIdx = blockBuff;
00389     }
00390   }
00394   void unsetAtomOverlp (void) {
00395 
00396     atoms           = NULL;
00397     cells           = NULL;
00398     overlapMap      = NULL;
00399     zone            = NULL;
00400     densityBlockIdx = NULL;
00401   }
00402 
00406   void writeInFileNew(std::ofstream& out) ;
00410   void writeSubDomainHeader(std::ofstream & out) ;
00414   void fillAtomInStructure(std::vector<std::vector<QCAtomIn> >& v) ;
00418    void updateInternalStructure() ;
00422    void printHeadOfInternalStructure();
00423 
00425 protected:
00426 
00427 
00429 private:
00430 
00431 
00435   void swap (int i, int j) {
00436     int tmp;
00437     tmp                           = overlapMap[i + LOCAL_IND_IDX];
00438     overlapMap[i + LOCAL_IND_IDX] = overlapMap[j + LOCAL_IND_IDX];
00439     overlapMap[j + LOCAL_IND_IDX] = tmp;
00440 
00441     tmp                            = overlapMap[i + REMOTE_IND_IDX];
00442     overlapMap[i + REMOTE_IND_IDX] = overlapMap[j + REMOTE_IND_IDX];
00443     overlapMap[j + REMOTE_IND_IDX] = tmp;
00444 
00445     tmp                             = overlapMap[i + REMOTE_ZONE_IDX];
00446     overlapMap[i + REMOTE_ZONE_IDX] = overlapMap[j + REMOTE_ZONE_IDX];
00447     overlapMap[j + REMOTE_ZONE_IDX] = tmp;
00448   }
00449 
00453     void sort (int offset) {
00454       bool sorted = false;
00455       while (!sorted) {
00456         sorted = true;
00457         for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]-1; ++i) {
00458           int cell = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00459           if (overlapMap[cell + REMOTE_IND_IDX] > 
00460               overlapMap[cell + IND_MAP_CELL + REMOTE_IND_IDX]) {
00461             sorted = false;
00462             swap (cell, cell + IND_MAP_CELL);
00463           }
00464         }
00465       }
00466     }
00467 
00470 public:
00471 protected:
00472 private:
00476   int domainId; 
00477   int globalDomainId; 
00481   const QCSystem * rootSystem;
00485   int nbAtomCoreShell1;
00489   int nbAOCoreShell1;
00493   int nbAtomShell2;
00497   int nbAOShell2;
00502   QCSubDomainZone * zone;
00507    int * globalIdx;
00511   int * densityBlockIdx;
00515   int densityBlockCount;
00519   int * overlapMap;
00523    int overlapMapSize;
00527   QCSymMatrix weight;
00533   QCFloat * bFactor;
00538   int   nbFilledOM;
00543   int   nbOccupiedOM;
00550   QCFloat * QCRestrict  orbitalOccupN;
00551 };
00552 
00556 extern ostream& operator << (ostream& out, const QCSubDomain& domain);
00557 
00558 
00559 #endif // QC_SUBDOMAIN_H_

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