QCSubDomain.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 <iomanip>
00021 #include <vector>
00022 
00023 #include "QCSubDomain.hpp"
00024 #include "QCTable.hpp"
00025 #include "QCMndoParam.hpp"
00026 #include "QCAm1Param.hpp"
00027 #include "QCPm3Param.hpp"
00028 
00029 
00030 
00034 QCSubDomain::QCSubDomain (const QCSystem *root) : QCSystem(),  QCModelMatrices(),
00035                                                   QCEnergies(),  domainId(-1),  globalDomainId(-1),rootSystem(root),
00036                                                   nbAtomCoreShell1(0),  nbAOCoreShell1(0),  nbAtomShell2(0),
00037                                                   nbAOShell2(0),  zone(NULL), globalIdx(NULL), densityBlockIdx(NULL),
00038                                                   densityBlockCount(0), overlapMap(NULL), overlapMapSize(0), weight(),
00039                                                   bFactor(NULL), nbFilledOM(0), nbOccupiedOM(0), orbitalOccupN(NULL)
00040 {}
00041 
00042 
00043 
00044 
00048 QCSubDomain::QCSubDomain (void) : QCSystem(), QCModelMatrices(), QCEnergies(), domainId(-1),globalDomainId(-1),
00049                                   rootSystem(NULL), nbAtomCoreShell1(0), nbAOCoreShell1(0),
00050                                   nbAtomShell2(0), nbAOShell2(0), zone(NULL), globalIdx(NULL),
00051                                   densityBlockIdx(NULL), densityBlockCount(0), overlapMap(NULL),
00052                                   overlapMapSize(0), weight(), bFactor(NULL), nbFilledOM(0),
00053                                   nbOccupiedOM(0), orbitalOccupN(NULL)
00054 {}
00055 
00056 
00057 
00058 
00062 QCSubDomain::~QCSubDomain (void) {
00063   QC_TRACE_END("BEGIN  QCSubDomain::~QCSubDomain");
00064   if (zone) {
00065     delete [] zone;  zone = NULL;
00066   }
00067   if (globalIdx) {
00068     delete [] globalIdx;  globalIdx = NULL;
00069   }
00070   if (densityBlockIdx) {
00071     delete [] densityBlockIdx;  densityBlockIdx = NULL;
00072   }
00073   if (overlapMap) {
00074     delete [] overlapMap;  overlapMap = NULL;
00075   }
00076   if (bFactor) {
00077     delete [] bFactor;   bFactor = NULL;
00078   }
00079   if (orbitalOccupN) {
00080     delete [] orbitalOccupN;  orbitalOccupN = NULL;
00081   }
00082   QC_TRACE_END("END    QCSubDomain::~QCSubDomain");
00083 }
00084 
00085 
00086 
00087 
00091 void
00092 QCSubDomain::allocStructures (int nbatoms, int nbao, int nbovrlp,  // nombre de sousdomaines
00093                               const int *sbdIds, const int *nbshared) {
00094   QC_TRACE_INIT("BEGIN  QCSubDomain::allocStructures for subdomain "<< domainId);
00095   QC_TRACE_INIT("              nbatoms  "<<nbatoms <<"  nbao  "<< nbao<<"  nbovrlp "<<nbovrlp);
00096  
00097   allocateAtoms(nbatoms);
00098   allocateCells(nbatoms);
00099   zone              = new QCSubDomainZone [nbatoms];
00100   globalIdx         = new int [nbatoms];
00101   densityBlockIdx   = new int [nbatoms];
00102   bFactor           = new QCFloat [nbao];
00103   orbitalOccupN     = new QCFloat [nbao];
00104   weight.setDimAndAllocate(nbatoms);
00105   //
00106   this->overlapMapSize = OVRLP_MAP_HEAD + nbovrlp * OVRLP_MAP_DECL;
00107   // temp
00108   int * offset = new int [nbovrlp];
00109 
00110   for (int i=0; i < nbovrlp; ++i) {
00111     offset[i]             = this->overlapMapSize;
00112     this->overlapMapSize += IND_MAP_HEAD + IND_MAP_CELL * nbshared[i];
00113   }
00114   
00115   this->overlapMap               = new int [this->overlapMapSize];
00116   this->overlapMap[NB_OVRLP_IDX] = nbovrlp;
00117 
00118   for (int i=0; i < nbovrlp; ++i) {
00119     
00120     this->overlapMap [OVRLP_MAP_HEAD + i            ] = offset[i];
00121     this->overlapMap [offset[i] + SBD_ID_IDX        ] = sbdIds[i];
00122     this->overlapMap [offset[i] + NB_SHARED_ATOM_IDX] = 0;
00123     this->overlapMap [offset[i] + NB_ATOM_SHELL1_IDX] = 0;
00124     this->overlapMap [offset[i] + NB_ATOM_SHELL2_IDX] = 0;
00125     this->overlapMap [offset[i] + NB_SHARED_AO_IDX  ] = 0;
00126     this->overlapMap [offset[i] + NB_AO_SHELL1_IDX  ] = 0;
00127     this->overlapMap [offset[i] + NB_AO_SHELL2_IDX  ] = 0;
00128   }
00129 
00130   delete [] offset;
00131   QC_TRACE_INIT("END  QCSubDomain::allocStructures");
00132   }
00133 
00139 void QCSubDomain::allocStructures (QCPartSubDomain & domain ){
00140  QC_TRACE_INIT("BEGIN QCSubDomain::allocStructures (QCPartSubDomain & domain )  ");
00141  //
00142  int nbAtoms = domain.getNbAtoms(), nbNeighbors= domain.getNbNeighbors(), nbAo =  domain.getLoad();
00143  //
00144  allocateAtoms(nbAtoms);
00145  allocateCells(nbAtoms);
00146  //
00147  zone              = new QCSubDomainZone [nbAtoms];
00148  globalIdx         = new int [nbAtoms];
00149  densityBlockIdx   = new int [nbAtoms];
00150  //
00151  weight.setDimAndAllocate(nbAtoms);
00152  //
00153  bFactor           = new QCFloat [nbAo];
00154  orbitalOccupN     = new QCFloat [nbAo];
00155  //
00156  this->overlapMapSize = OVRLP_MAP_HEAD + nbNeighbors * OVRLP_MAP_DECL;
00157  int * offset  = new int [nbNeighbors];
00158  //
00159  for(int i = 0 ; i < domain.getNbNeighbors(); ++i){
00160    offset[i]             = this->overlapMapSize;
00161    this->overlapMapSize += IND_MAP_HEAD + IND_MAP_CELL * ( domain.getNbAtomsShell1(i) + domain.getNbAtomsShell2(i));
00162  }
00163  this->overlapMap               = new int [this->overlapMapSize];
00164  this->overlapMap[NB_OVRLP_IDX] = nbNeighbors;
00165  
00166  for (int i=0; i < nbNeighbors; ++i) {
00167    this->overlapMap [OVRLP_MAP_HEAD + i            ] = offset[i];
00168    this->overlapMap [offset[i] + SBD_ID_IDX        ] = domain.getNeighborIdx(i) ;
00169    this->overlapMap [offset[i] + NB_SHARED_ATOM_IDX] = 0 ;
00170    this->overlapMap [offset[i] + NB_ATOM_SHELL1_IDX] = domain.getNbAtomsShell1(i);
00171    this->overlapMap [offset[i] + NB_ATOM_SHELL2_IDX] = domain.getNbAtomsShell2(i);
00172    this->overlapMap [offset[i] + NB_SHARED_AO_IDX  ] = 0;
00173    this->overlapMap [offset[i] + NB_AO_SHELL1_IDX  ] = 0;
00174    this->overlapMap [offset[i] + NB_AO_SHELL2_IDX  ] = 0;
00175  }
00176 
00177  delete [] offset;
00178  //
00179 #ifdef DEBUG_INIT_SUBDOMAIN
00180  this->printHeadOfInternalStructure();
00181 #endif
00182  QC_TRACE_INIT("END   QCSubDomain::allocStructures (QCPartSubDomain & domain )  ");
00183 }
00184 
00188 void
00189 QCSubDomain::addOverlpIndirectionInfo (int remoteSbdId,    int atomIdx,
00190                                        int remoteAtomIdx,   QCSubDomainZone remoteZone) {
00191  
00192  
00193  int offset = this->getMapOffset(remoteSbdId);
00194 
00195 #ifndef QC_NO_DEBUG
00196   assert(offset >= 0);
00197 #endif
00198   
00199     
00200   int& nbCells   = overlapMap[offset + NB_SHARED_ATOM_IDX];
00201   int cellOffset = ( offset + IND_MAP_HEAD +  nbCells * IND_MAP_CELL);
00202   //  
00203   overlapMap [cellOffset + LOCAL_IND_IDX  ] = atomIdx;
00204   overlapMap [cellOffset + REMOTE_IND_IDX ] = remoteAtomIdx;
00205   overlapMap [cellOffset + REMOTE_ZONE_IDX] = remoteZone;
00206 
00207   nbCells++;
00208 } 
00209 
00213 template <class TPParam>
00214 void 
00215 QCSubDomain::setSystemNumbers (const TPParam *QCRestrict parameters) {
00216   QC_TRACE_INIT("BEGIN QCSubDomain::setSystemNumbers  ");
00217   //
00218   nbAtomicOrbitals = densityBlockCount = nbElectrons = 0 ; 
00219   nbAtomShell2 = nbAtomCoreShell1 = 0 ;
00220   //
00221   for (int i=0; i < nbAtoms; i++) {
00222     
00227     setFirstAOAt(i, nbAtomicOrbitals);
00228     setDensityBlockIdx(i, densityBlockCount);
00229     int nbAO = parameters[getParamIndexAt(i)].getNbAO();
00230     nbAtomicOrbitals  += nbAO;
00231     densityBlockCount += (nbAO * (nbAO + 1)) / 2;
00232     nbElectrons       += parameters[getParamIndexAt(i)].getNbValenceElec();
00233 //      if (zone[i] == QC_SHELL2) {
00234 //       ++nbAtomShell2;
00235 //       nbAOShell2 += parameters[getParamIndexAt(i)].getNbAO(); 
00236 //     } 
00237 //     else {
00238 //       ++nbAtomCoreShell1;
00239 //       nbAOCoreShell1  += parameters[getParamIndexAt(i)].getNbAO();
00240 //     }
00241   }
00245   nbElectrons      -= systemCharge * (nbElectrons / rootSystem->getNbElectrons());
00246   spinMultiplicity  = rootSystem->getSpinMultiplicity();
00247 #if defined TRACE_OUT
00248   this->printHeadOfInternalStructure();
00249 #endif
00250   QC_TRACE_INIT("END   QCSubDomain::setSystemNumbers  ");
00251 }
00252 
00256 void QCSubDomain::fillWeight (const QCGeneralData& data) {
00257 //
00258   int k, i, j;
00259   int localIdxA, localIdxB;
00260   
00261   for (i = 0; i < getNbAtoms(); ++i) {
00262     for (j = 0; j <= i; ++j) {
00263       weight[i][j] = 1;
00264     }
00265   }
00266 
00267   // On parcours la liste des SD recouvrants
00268   for (k = 0; k < getNbOverlpSD(); ++k) {
00269     int offset = overlapMap[OVRLP_MAP_HEAD + k];
00270     
00271     //On parcours les atomes de la zone de recouvrement
00272     for (i = 0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00273       int cellA = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00274       localIdxA = overlapMap[cellA + LOCAL_IND_IDX]; 
00275 
00276       // On selectionne les couples d atomes qui appartiennent a la zone
00277       // de recouvrement
00278       for (j = 0; j <= i; ++j) {
00279         int cellB = offset + IND_MAP_HEAD + j * IND_MAP_CELL;
00280         localIdxB = overlapMap[cellB + LOCAL_IND_IDX]; 
00281 
00282         //Attention, la matrice est symetrique.
00283         if (localIdxA > localIdxB) {
00284           weight[localIdxA][localIdxB] += 1;
00285         
00286         } else {
00287           weight[localIdxB][localIdxA] += 1;
00288         }
00289       }
00290     }    
00291   }
00292   
00293   // Et maintenant, on inverse.
00294   for (i = 0; i < getNbAtoms(); ++i) {
00295     for (j = 0; j <= i; ++j) {
00296       // Le coef de partitionnement est nul des que
00297       // l'on est dans la zone de cut
00298       if (data.getPartitionType() == QC_DIXON_PART) {
00299         
00300         if (zone[i] == QC_SHELL2 || zone[j] == QC_SHELL2) {
00301           weight[i][j] = QC_ZERO;         
00302       
00303         } else if (zone[i] == QC_SHELL1 && zone[j] == QC_SHELL1) {
00304           weight[i][j] = QC_ZERO;
00305         
00306         } else {
00307           weight[i][j] = QC_ONE / weight[i][j];
00308         }
00309       
00310       } else if (data.getPartitionType() == QC_STANDARD_PART) {
00311         
00312         if (zone[i] == QC_SHELL2 || zone[j] == QC_SHELL2) {
00313           weight[i][j] = QC_ZERO;
00314         
00315         } else {
00316           weight[i][j] = QC_ONE / weight[i][j];
00317         }      
00318       
00319       }
00320     }
00321   }
00322 }
00323 
00327 template <class TPParam>
00328 void
00329 QCSubDomain::fillDensityArray (const TPParam * QCRestrict parameters,
00330                                QCFloat  * QCRestrict &    densityArray,
00331                                bool                       alloc) {
00332   
00333   // Le tableau de parametres
00334   const TPParam * QCRestrict param;
00335 
00336   // La matrice densite
00337   const QCSymMatrix& density = getDensityP().getMatrix();
00338 
00339   // les indices
00340   int i, mu, nu, firstAO, lastAO, blockIdx; 
00341 
00342   if (alloc) {
00343     densityArray = new QCFloat [densityBlockCount];
00344   }
00345 
00346   blockIdx = 0;
00347   
00348   for (i=0; i < nbAtoms; ++i) {
00349     param   = &parameters[getParamIndexAt(i)];
00350     firstAO = getFirstAOAt(i);
00351     lastAO  = firstAO + param->getNbAO();
00352 
00353     for (mu = firstAO; mu < lastAO; ++mu) {
00354       for (nu = firstAO; nu <= mu; ++nu) {
00355         densityArray [blockIdx] = density [mu][nu];
00356         ++blockIdx;
00357       }
00358     }
00359   }
00360 }
00361 
00367 template <class TPParam>
00368 void
00369 QCSubDomain::initWeightedDensity (const TPParam * QCRestrict parameters) {
00370 
00371   // La matrice densite
00372    QCSymMatrix& QCRestrict densityP = getDensityP().getMatrix();
00373 
00374   // La matrice densit�pender�
00375   QCSymMatrix& QCRestrict weightedDensityP = 
00376     getWeightedDensityP()->getMatrix();
00377 
00378   // Les parametres
00379   const TPParam * QCRestrict paramA;
00380   const TPParam * QCRestrict paramB;
00381 
00382   // Les orbitales atomiques
00383   int firstAOofA, firstAOofB;
00384   int lastAOofA, lastAOofB;
00385 
00386   // Les poids
00387   QCFloat weightAA, weightAB;
00388 
00389   // les indices
00390   int mu, nu, lda;
00391   int a, b;
00392 
00393   
00394   for (a = 0; a < getNbAtoms(); ++a) {
00395     
00396     paramA     = &parameters[getParamIndexAt(a)];
00397     firstAOofA = getFirstAOAt(a);
00398     lastAOofA  = firstAOofA + paramA->getNbAO();    
00399     
00400     for (b = 0; b < a; ++b) {
00401       
00402       paramB     = &parameters[getParamIndexAt(b)];
00403       firstAOofB = getFirstAOAt(b);
00404       lastAOofB  = firstAOofB + paramB->getNbAO();    
00405 
00406       weightAB         = weight[a][b];
00407 
00408       for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00409         for (lda = firstAOofB; lda < lastAOofB; ++lda) {
00410           weightedDensityP[mu][lda] = weightAB * densityP[mu][lda];
00411         }
00412       }
00413     }
00414     
00415     // On n'oublie pas le cas ou b == a.
00416     weightAA = weight[a][a];
00417     
00418     for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00419       for (nu = firstAOofA; nu <= mu; ++nu) {
00420         weightedDensityP[mu][nu] =  weightAA * densityP[mu][nu];
00421       }
00422     }
00423   }
00424 }
00425 
00429 template <class TPParam>
00430 void
00431 QCSubDomain::computeBFactor (const TPParam * QCRestrict parameters) {
00432 
00433   // La matrice des vecteurs propres
00434   QCMatrix& QCRestrict localEigenVectCt = getEigenVectCt().getMatrix();
00435 
00436   // Les premiers index d'orbitales
00437   int firsAOofB, lastAOofB;
00438 
00439   // Le parametre courant
00440   const TPParam * QCRestrict paramB;
00441 
00442   // Le poid de l'atome courant
00443   QCFloat weightBB;
00444 
00445   // Les indices
00446   int b, i, mu;
00447 
00448 
00449   for (i = 0; i < getNbAtomicOrbitals(); ++i) {
00450     bFactor[i] = QC_ZERO;
00451     for (b = 0; b < getNbAtoms(); ++b) {
00452       weightBB  = weight[b][b];
00453       paramB    = &parameters[getParamIndexAt(b)];
00454       firsAOofB = getFirstAOAt(b);
00455       lastAOofB = firsAOofB + paramB->getNbAO();
00456       
00457       for (mu = firsAOofB; mu < lastAOofB; ++mu) {
00466         bFactor[i] += weightBB * FQCPow<2>(localEigenVectCt[i][mu]);
00467       }
00468     }
00469     // cout << "bfactor[" << i << "] = " << bFactor[i] << endl;
00470   }
00471 }
00472 
00473 //                                                                            //
00474 //----------------------------------------------------------------------------//
00475 // Private Method: traceProductCoreAndOvlp
00476 // 
00477 // Description: Trace du produit de 2 matrices uniquement sur les zones core et Buffer1.
00478 //
00479 template <class TPParam>
00480 QCFloat 
00481 QCSubDomain::traceProductCoreShell1 (const TPParam * QCRestrict parameters,
00482                                      QCSymMatrix&    QCRestrict fockMatrix,
00483                                      QCSymMatrix&    QCRestrict densityMatrix) {
00484 
00485 
00486   // leurs parametres
00487   const TPParam * QCRestrict parameterA;
00488   const TPParam * QCRestrict parameterB;
00489 
00490 
00491   // Les premiers index d'orbitales de A et B dans la matrice complete du
00492   // systeme et l'index du premier successeur de A tjs dans la 
00493   // matrice complete.
00494   int firstAOofA, firstAOofB, lastAOofA, lastAOofB;
00495 
00496   // valeurs intermediaire
00497   register QCFloat interVal, interVal2, interVal3;
00498   
00499   //
00500   int i, j, mu, nu, lda;
00501 
00502   // La trace a retourner
00503   QCFloat traceFP = QC_ZERO;
00504 
00505   for (i = 0; i < getNbAtoms(); ++i) {
00506     
00507     if (getZone(i) != QC_SHELL2) {
00508       parameterA = &parameters[getParamIndexAt(i)];
00509       firstAOofA = getFirstAOAt(i);
00510       lastAOofA  = firstAOofA + parameterA->getNbAO();
00511 
00512       for (j = 0; j < i; ++j) {
00513 
00514         if (getZone(j) == QC_SHELL2) {
00515           continue ;
00516         }
00517           //    if (getZone(j) != QC_SHELL2) {
00518         parameterB = &parameters[getParamIndexAt(j)];
00519         firstAOofB = getFirstAOAt(j);
00520         lastAOofB  = firstAOofB + parameterB->getNbAO();
00521 
00522         // Si on est dans (overlap inter overlap) ou dans 
00523         // (overlap inter core) ou dans (core inter core).
00524         interVal = QC_ZERO;
00525         for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00526           for (lda = firstAOofB; lda < lastAOofB; ++lda) {
00527             //
00528             // Attention, E = 1/2 (trace(H+F) D).
00529             // Le coef 1/2 n est pas utilise car ici, on traite l interaction A-B.
00530             // Le traitement symetrique de B-A donnerait la meme valeur.
00531             // Donc, on omet le facteur 1/2 et on n a pas besoin de traiter B-A.
00532             interVal += (fockMatrix[mu][lda] * densityMatrix[mu][lda]);
00533           }
00534         }
00535 
00536         // Si on est dans (overlap inter core), on doit virer completement
00537         // la contribution a la trace des atomes parcourus.
00538         if (getZone(i) == QC_CORE || getZone(j) == QC_CORE) {
00539           traceFP += interVal;
00540         }  
00541         // Si on est dans (overlap inter overlap), on multiplie
00542         //       par le poids avant de virer.
00543         else {
00544           // On est dans (overlap inter overlap)
00545           traceFP += weight[i][j] * interVal;
00546         }
00547       }
00548       //    }
00549       // Traitement des termes diagonaux.
00550       interVal  = QC_ZERO;
00551       interVal2 = QC_ZERO;
00552       for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00553         for (nu = firstAOofA; nu < mu; ++nu) {
00554           interVal += (fockMatrix[mu][nu] * densityMatrix[mu][nu]);
00555         }
00556         interVal2 += (fockMatrix[mu][mu] * densityMatrix[mu][mu]);
00557       }
00558       //
00559       // On ne compte les termes diagonaux qu'une seule fois car ils
00560       // n'ont pas de symetrique.
00561       interVal3 = interVal + QC_HALF * interVal2;
00562       
00563       // Si l'atome est dans la zone de core, on doit virer completement
00564       // sa contribution a la trace.
00565       if (getZone(i) == QC_CORE) {
00566         traceFP += interVal3;
00567       }
00568       // S'il apartient a la zone d overlap, on mutiplie par 1/n 
00569       //     avant de virer sa contribution a la trace.
00570       else {
00571         //  getZone(i) == QC_SHELL1
00572         traceFP += weight[i][i] * interVal3;
00573       }
00574     }
00575   }
00576   return traceFP;
00577 }
00578 
00582 void QCSubDomain::modifyMatricesAndEnergy (const QCFloat& lambdaOpt) {
00583 
00584   QCFloat energyFockTild;
00585   QCFloat energyTild1e, c1;
00586   
00587   c1 = QC_ONE - lambdaOpt ; 
00588   fockFTild->getMatrix().scale(c1);
00589   fockFTild->getMatrix().axpy(lambdaOpt, fockF.getMatrix());
00590 
00591   densityPTild->getMatrix().scale(c1);
00592   densityPTild->getMatrix().axpy(lambdaOpt, densityP.getMatrix());
00593 
00594   // On remplace F par F_tilde
00595   fockF.copyMatrix(*fockFTild);
00596   //
00597   // Calcul des energies
00598   //
00599   energyFockTild = ( QC_HALF * 
00600                      fockFTild->getMatrix().traceProduct(densityPTild->getMatrix()) );
00601   energyTild1e   = ( QC_HALF * 
00602                      hamiltonH.getMatrix().traceProduct(densityPTild->getMatrix()) );
00603   //
00604   elecEnergyETild   = energyTild1e + energyFockTild;
00605   elecEnergyETild1e = energyTild1e;
00606 }
00607 
00611 ostream& operator << (ostream& out, const QCSubDomain& domain) {
00612 
00613   out << " * Atoms:" << endl;    
00614   out << "{" << endl << " LOCAL/GLOBAL:" << " (           X" 
00615       << ",            Y,            Z), TYPE,    ZONE, (SD i: LOCAL)" 
00616       << endl << endl;
00617 
00618   for (int i=0; i < domain.getNbAtoms(); ++i) {
00619     QCPoint3D current;
00620     
00621     domain.getPointAt(i, current);
00622 
00623     out << " " << setw(5) << i << "/ " 
00624         << setw(5) << domain.getGlobalIdx(i) << ": "
00625         << current                << ",    "
00626         << table[domain.getType(i)]          << ",  "
00627         << QCCommon::domainZones.find(domain.getZone(i))->second
00628         << ", ";
00629     
00630     for (int j=0; j < domain.overlapMap[NB_OVRLP_IDX]; ++j) {
00631       int offset    = domain.overlapMap[OVRLP_MAP_HEAD + j];
00632       int domainId  = domain.overlapMap[offset + SBD_ID_IDX]; 
00633       int remoteIdx = domain.getRemoteIdx(domainId, i);
00634       
00635       if (remoteIdx >= 0) {
00636         out << "(sd" << setw(2) << domainId  << ": " << setw(5) << remoteIdx  << ") ";
00637       }
00638     }
00639 
00640     out << endl;
00641   }
00642   out << "}" << endl;
00643   out << " * Overlapment:" << endl;
00644   out << "{" << endl;
00645   
00646   for (int j=0; j < domain.overlapMap[NB_OVRLP_IDX]; ++j) {
00647     int offset = domain.overlapMap[OVRLP_MAP_HEAD + j];
00648 
00649     out << "* domain " << setw(2) << domain.getId() 
00650         << " <-> domain " << setw(2) << domain.overlapMap[offset + SBD_ID_IDX]
00651         << "  (" << setw(4) << domain.overlapMap[offset + NB_SHARED_ATOM_IDX]
00652         << " atoms):" << endl
00653         << "   {" << endl;
00654     
00655     out << "  (LOCAL, GlOBAL  ZONE) -> (REMOTE,  RZONE)" << endl << endl;
00656     for (int i=0; i < domain.overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00657       
00658       int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00659       int localIdx   = domain.overlapMap[cellOffset + LOCAL_IND_IDX];
00660       int remoteIdx  = domain.overlapMap[cellOffset + REMOTE_IND_IDX];
00661       int remoteZone = domain.overlapMap[cellOffset + REMOTE_ZONE_IDX];
00662 
00663       out << "  (" << setw(5)  << localIdx << "   "<< domain.globalIdx[localIdx]
00664           << ", "  << QCCommon::domainZones.find(domain.getZone(localIdx))->second
00665           << ") -> ( "  << setw(5)  << remoteIdx          << ", " 
00666           << QCCommon::domainZones.find(static_cast<QCSubDomainZone>(remoteZone))->second
00667           << ")";
00668       out << endl;
00669     }
00670     out << "   }" << endl << endl;
00671   }
00672   out << "}" << endl;
00673 
00674   out << "\n * Nb Atoms:          " << domain.getNbAtoms()
00675       << " (core+shell1: "          << domain.getNbAtomCoreShell1()
00676       << ",   shell2: "             << domain.getNbAtomShell2()
00677       << ")"
00678       << "\n * System Charge:     " << domain.getSystemCharge()
00679       << "\n * Spin Multiplicity: " << domain.getSpinMultiplicity()
00680       << "\n * NbAtomicOrbitals:  " << domain.getNbAtomicOrbitals()
00681       << " (core+shell1: "          << domain.getNbAOCoreShell1()
00682       << ",   shell2: "             << domain.getNbAOShell2()
00683       << ")"
00684       << "\n * NbElectrons:       " << domain.getNbElectrons() 
00685       << endl;
00686   
00687   return out;
00688 }
00689 
00690 void QCSubDomain::writeInFileNew(std::ofstream& out){
00691 
00692   QC_TRACE_OUT("BEGIN QCSubDomain::writeInFileNew");
00693   //
00694 #if defined TRACE_OUT
00695   this->printHeadOfInternalStructure();
00696 #endif
00697   int nbSD, numCore,  numOverlap, numCut, n, j;
00698   //
00699   nbSD = this->overlapMap[NB_OVRLP_IDX];
00700   numCore =  numOverlap = numCut = 0 ;
00702   for (  n = 0 ; n < this->getNbAtoms() ; ++n){
00703     if (this->zone[n] == QC_CORE){
00704       ++numCore ; }
00705     else if (this->zone[n] == QC_SHELL1){
00706       ++numOverlap ; }
00707     else {
00708       ++numCut ; }   
00709   }
00710   out << "Sub Domain number : " << this->domainId <<std::endl ;
00711   out << " "<< nbSD <<"    " << numCore <<"  "<< numOverlap << "  " << numCut<<std::endl ;
00712   out << "   CORE ZONE     : " <<std::endl ;
00713   if (numCore !=0){
00714     out <<  "      ";
00715     for (  n = 0 ; n < this->getNbAtoms()  ; ++n){
00716       if (this->zone[n] == QC_CORE){
00717         out <<  this->globalIdx[n] << " " ;
00718       }
00719     }
00720     out << std::endl ;
00721   }
00722 #ifdef DEBUG_INIT_SUBDOMAIN
00723   this->printHeadOfInternalStructure();
00724 #endif
00725   int offset , nbShared, nAt ,localIdx;
00726   out << "   OVERLAP ZONE  : " <<std::endl ;
00727   for (  n = 0 ; n < nbSD ; ++n){
00728     offset   = this->overlapMap[OVRLP_MAP_HEAD + n*OVRLP_MAP_DECL] ;
00729     nbShared = this->overlapMap[offset + NB_SHARED_ATOM_IDX] ;
00730     nAt      = this->overlapMap[offset + NB_ATOM_SHELL1_IDX] ;
00731     if (nAt == 0 ) { continue;}
00732     out << "      " << this->overlapMap[offset + SBD_ID_IDX] <<"  "<< nAt << "      ";
00733     //
00734     offset += IND_MAP_HEAD + LOCAL_IND_IDX;
00735     for( j = 0 ; j <  nbShared ; ++j){
00736       localIdx = overlapMap[offset + j*IND_MAP_CELL];
00737       if( this->zone[localIdx] == QC_SHELL1){
00738         out << this->globalIdx[localIdx]<< " " ;
00739       }
00740     }
00741     out << std::endl;
00742   }
00743   //  std::cout <<"   CUT ZONE      : "  <<std::endl ;
00744   out << "   CUT ZONE      : "  <<std::endl ;
00745   for ( n = 0 ; n < nbSD ; ++n){
00746     offset   = this->overlapMap[OVRLP_MAP_HEAD  + n*OVRLP_MAP_DECL] ;
00747     nbShared = this->overlapMap[offset + NB_SHARED_ATOM_IDX] ;
00748     nAt      = this->overlapMap[offset + NB_ATOM_SHELL2_IDX] ;
00749     //    std::cout << "    SBD_ID_IDX: "<< this->overlapMap[offset + SBD_ID_IDX] <<"  nAt "<< nAt << std::endl;
00750     if (nAt == 0 ) { continue;}
00751     //    std::cout <<"   ....write : " << std::endl;
00752     out << "      " << this->overlapMap[offset + SBD_ID_IDX] <<"  "<< nAt << "      ";
00753     offset += IND_MAP_HEAD + LOCAL_IND_IDX ;
00754     for( j = 0; j <  nbShared ; ++j){
00755       localIdx = overlapMap[offset + j*IND_MAP_CELL];
00756 //       std::cout <<"  localIdx : " <<localIdx << " globalIdx[localIdx]  " <<  this->globalIdx[localIdx]
00757 //              << "  zone[localIdx] " <<this->zone[localIdx]<< std::endl;
00758       if( this->zone[localIdx] == QC_SHELL2){
00759         out << this->globalIdx[localIdx]<< " " ;
00760       }
00761     }
00762     out << std::endl;
00763   }
00764   //
00765   QC_TRACE_OUT("END   QCSubDomain::writeInFileNew");
00766 }
00767 
00769  void QCSubDomain::updateInternalStructure(){
00771    int n,i, offset, localIdx,cellOffset, nbSHELL1, nbSHELL2, nbOASHELL1, nbOASHELL2;
00773    QC_TRACE_INIT("BEGIN QCSubDomain::updateInternalStructure "<< this->domainId << " GlobalIdx " <<globalDomainId);
00774 #ifdef DEBUG_INIT_SUBDOMAIN
00775    std::cout<< std::endl<< std::endl<< " HEAD updated "<< std::endl
00776             <<"  Number Of SDs : "<<this->overlapMap[NB_OVRLP_IDX]<< std::endl;
00777    this->printHeadOfInternalStructure();
00778 #endif
00779    // 
00780    //
00781    nbOASHELL1  = nbOASHELL2 = 0 ;
00782    int oa;
00784    //
00785    for ( n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00786      
00787      offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00789      nbSHELL1 = nbSHELL2 = 0 ; nbOASHELL1  = nbOASHELL2 = 0 ;
00790      for ( i=0; i < this->overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00791        
00792        cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00793        localIdx   = this->overlapMap[cellOffset + LOCAL_IND_IDX];
00794        //       oa         = (table[this->QCSystem::getType(localIdx)] == "H" )? 1 :4 ;
00795        oa = (this->QCSystem::getType(localIdx) == 1 /* H*/)? 1 : 4 ; 
00796        if (this->zone[localIdx] == QC_SHELL1){
00797          ++nbSHELL1 ;  nbOASHELL1 += oa ;
00798        }
00799        else if (this->zone[localIdx] == QC_SHELL2){
00800          ++nbSHELL2; nbOASHELL2 += oa ;
00801        }   
00802      }
00803      this->overlapMap[offset + NB_ATOM_SHELL1_IDX] = nbSHELL1 ;
00804      this->overlapMap[offset + NB_AO_SHELL1_IDX]   = nbOASHELL1 ;
00805      this->overlapMap[offset + NB_ATOM_SHELL2_IDX] = nbSHELL2 ;
00806      this->overlapMap[offset + NB_AO_SHELL2_IDX]   = nbOASHELL2 ;
00807      this->overlapMap[offset + NB_SHARED_ATOM_IDX] = nbSHELL1+ nbSHELL2;
00808      this->overlapMap[offset + NB_SHARED_AO_IDX]   = nbOASHELL1+ nbOASHELL2;
00809    }
00810    //     oa = (table[this->QCSystem::getType(at)] == "H" )? 1 : 4 ; 
00811 
00812    //
00813    this->nbAtomCoreShell1 = this->nbAOCoreShell1 = this->nbAtomShell2 = this->nbAOShell2 = 0 ;   
00814    for (int at =0; at < this->QCSystem::getNbAtoms() ;++at){
00815      //     oa = (table[this->QCSystem::getType(at)] == "H" )? 1 : 4 ; 
00816      oa = (this->QCSystem::getType(at) == 1 )? 1 : 4 ; 
00817      if (this->zone[at] == QC_SHELL2){
00818        ++this->nbAtomShell2 ;  this->nbAOShell2  += oa;
00819      }
00820      else {
00821        this->nbAOCoreShell1 += oa; ++this->nbAtomCoreShell1;
00822      }
00823    }
00824 #ifdef DEBUG_INIT_SUBDOMAIN
00825    this->printHeadOfInternalStructure();
00826 #endif
00827 
00828    QC_TRACE_INIT("END  QCSubDomain::updateInternalStructure");
00829  }
00830 void QCSubDomain::printHeadOfInternalStructure(){
00832  int n, offset;
00836  QC_TRACE("BEGIN QCSubDomain::printHeadOfInternalStructure  " << this->overlapMap[NB_OVRLP_IDX] );
00837  std::cout <<" ############################################################################################" <<std::endl;
00838  std::cout <<" Head Of internal structure of subdomain "<< domainId << " numGlob : " << globalDomainId <<std::endl
00839            << " nbAtomCoreShell1 "<< this->nbAtomCoreShell1 <<std::endl
00840            << " nbAOCoreShell1   "<< this->nbAOCoreShell1   <<std::endl
00841            << " nbAtomShell2     "<< this->nbAtomShell2     <<std::endl
00842            << " nbAOShell2       "<< this->nbAOShell2       <<std::endl
00843            <<std::endl ;
00844  for ( n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00845    offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00846    std::cout << "SD num "<< n << " offset : " <<offset << " HEAD : "  << std::endl;
00847    std::cout<<  "  SBD_ID_IDX : " << this->overlapMap[offset + SBD_ID_IDX ] << "  "  << std::endl;
00848    std::cout<<  "       NB_SHARED_ATOM : " << this->overlapMap[offset + NB_SHARED_ATOM_IDX ] << "  "  << std::endl;
00849    std::cout<<  "       NB_ATOM_SHELL1 : " << this->overlapMap[offset + NB_ATOM_SHELL1_IDX ] << "  "  << std::endl;
00850    std::cout<<  "       NB_ATOM_SHELL2 : " << this->overlapMap[offset + NB_ATOM_SHELL2_IDX ] << "  "  << std::endl;
00851    std::cout<<  "       NB_SHARED_AO   : " << this->overlapMap[offset + NB_SHARED_AO_IDX   ] << "  "  << std::endl;
00852    std::cout<<  "       NB_AO_SHELL1   : " << this->overlapMap[offset + NB_AO_SHELL1_IDX   ] << "  "  << std::endl;
00853    std::cout<<  "       NB_AO_SHELL2   : " << this->overlapMap[offset + NB_AO_SHELL2_IDX   ] << "  "  << std::endl;
00854    std::cout << std::endl;
00855  }
00856  std::cout <<" ############################################################################################" <<std::endl;
00857  QC_TRACE("END   QCSubDomain::printHeadOfInternalStructure  ");
00858 }
00862 void QCSubDomain::writeSubDomainHeader(std::ofstream & out){
00863   QC_TRACE("BEGIN QCSubDomain::writeSubDomainHeader  ");
00864   int offset =0 ;
00865   int nbAtom , nbAO ;
00866   //
00867   nbAtom = this->nbAtomCoreShell1 + this->nbAtomShell2  ;
00868   nbAO   = this->nbAOCoreShell1   + this->nbAOShell2    ;
00869   out << "    "<<  this->domainId << "      " <<  nbAtom <<"      "  << nbAO << "      " 
00870       << this->overlapMap[NB_OVRLP_IDX]     << "      " ;
00871   //  std::cout << "    "<<  this->domainId << "      " <<  nbAtom <<"      "  << nbAO << "      " 
00872   //        << this->overlapMap[NB_OVRLP_IDX]     << "      " <<std::endl;
00873   for ( int n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00874     offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00875     out << this->overlapMap[offset + SBD_ID_IDX] <<"      "<< this->overlapMap[offset + NB_SHARED_ATOM_IDX] << "      ";
00876   }
00877   out << std::endl;
00878   QC_TRACE("END   QCSubDomain::writeSubDomainHeader  ");
00879 }
00880 void QCSubDomain::fillAtomInStructure(std::vector<std::vector<QCAtomIn> >& atomInZone){
00881   //
00882   QC_TRACE("BEGIN QCSubDomain::fillAtomInStructure  ");
00883   int n ; 
00884   QCAtomIn current;
00885   //
00886   current.numDomain = this->domainId ;
00887   for (  n = 0 ; n < this->getNbAtoms()  ; ++n){
00888     current.localNum = n;
00889     current.typeZone = this->zone[n];
00890     //
00891     (atomInZone[this->globalIdx[n]]).push_back(current);
00892   }
00893   QC_TRACE("END   QCSubDomain::fillAtomInStructure  ");
00894 }
00895 
00899 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::setSystemNumbers);
00900 
00901 QCPARAMETER_METH_EXPL_INST_PARAM(void QCSubDomain::fillDensityArray,
00902                                  TWO_PARAMS(QCFloat * QCRestrict &, bool));
00903 
00904 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::initWeightedDensity);
00905 
00906 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::computeBFactor);
00907 
00908 QCPARAMETER_METH_EXPL_INST_PARAM(QCFloat QCSubDomain::traceProductCoreShell1,
00909                                  TWO_PARAMS(QCSymMatrix& QCRestrict, QCSymMatrix& QCRestrict));

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