QCMDSystem.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 <sstream>
00021 #include <iomanip>
00022 #include <algorithm>
00023 #include <float.h>
00024 
00025 #include "QCMacro.hpp"
00026 #include "QCCommon.hpp"
00027 #include "QCReader.hpp"
00028 #include "QCMndo.hpp"
00029 #include "QCAm1.hpp"
00030 #include "QCPm3.hpp"
00031 #include "QCMndoParam.hpp"
00032 #include "QCAm1Param.hpp"
00033 #include "QCPm3Param.hpp"
00034 #include "QCMDSystem.hpp"
00035 #include "QCDistMDSystem.hpp"
00036 #include "QCSCFAlgo.hpp" 
00037 #include "QCDCAlgo.hpp"
00038 #include "QCManager.hpp"
00039 #include "QCIntgReader.hpp"
00040 #include "QCLinearPartitioner.hpp"
00041 #include "QCKimikaPartitioner.hpp"
00042 #include "QCTools.hpp"
00043 
00044 #if defined (QC_TRACE_IN_FILE)  
00045 #include "QCTrace.hpp"
00046 #endif
00047 
00050 template <class TPSolver>
00051 QCMDSystem<TPSolver>::QCMDSystem (void) : QCSystem(), QCEnergies(),  totalNbDomains(0), 
00052                                           nbDomains(0), domains(NULL),offset(0),  
00053                                           fermiFCall(true),_partitionType(QC_STANDARD_PART) 
00054 {}
00055 
00059 template <class TPSolver>
00060 QCMDSystem<TPSolver>::~QCMDSystem (void) {
00061   QC_TRACE_END("BEGIN QCMDSystem<TPSolver>::~QCMDSystem");
00062   if (domains != NULL) {
00063     delete [] domains;     domains = NULL;
00064   }
00065   QC_TRACE_END("END   QCMDSystem<TPSolver>::~QCMDSystem");
00066 }
00067 
00068 template <class TPSolver>
00069 bool QCMDSystem<TPSolver>::checkData(const QCGeneralData& data){
00070   QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00071 
00072   if(!data.isDCComputation()){
00073     std::cerr << "   Wrong algorithm for the SCF computation must have the DC_ prefix" << std::endl ;
00074     std::cerr << "       The value is "<< data.getMainAlgorithmStr() <<std::endl;
00075     exit(EXIT_FAILURE);
00076   }
00077    if(data.isReadPartitionFromFile()){
00078      // toDO
00079      std::cerr << "  Check todo in QCMDSystem<TPSolver>::checkData" << std::endl ;
00080    }
00081   QC_TRACE_INIT("END   QCMDSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00082    return true ;
00083 }
00087 template <class TPSolver>
00088 void 
00089 QCMDSystem<TPSolver>::addAtom (const QCPoint3D& coords, int type, int globalIdx,
00090                                int domainIdx, QCSubDomainZone zone, int atomIdx) {
00091   
00092   domains[domainIdx].setPointAt(atomIdx, coords);
00093   domains[domainIdx].setType(atomIdx, type);
00094   domains[domainIdx].setZone(atomIdx, zone);
00095   domains[domainIdx].setGlobalIdx(atomIdx, globalIdx);
00096   domains[domainIdx].setParamIndexAt(atomIdx, 
00097                                      getRootSystem().getParamIndexAt(globalIdx));
00098 }
00099 
00103 template <class TPSolver>
00104 void
00105 QCMDSystem<TPSolver>::fillMaps (int  nbsbd, const int * domainIdxs, const QCSubDomainZone * zones,
00106                                 const int * atomIdxs) {
00107   if (nbsbd > 1) {
00108     for (int i=0; i < nbsbd-1; ++i) {
00109       for (int j=i+1; j < nbsbd; ++j) {
00110         if (zones[i] != QC_CORE && zones[j] != QC_SHELL2) {
00111           domains[domainIdxs[i]].addOverlpIndirectionInfo(domainIdxs[j], atomIdxs[i],
00112                                                           atomIdxs[j],zones[j]);
00113         }
00114         if (zones[j] != QC_CORE && zones[i] != QC_SHELL2) {
00115           domains[domainIdxs[j]].addOverlpIndirectionInfo(domainIdxs[i], atomIdxs[j],
00116                                                           atomIdxs[i],zones[i]);
00117           
00118         }
00119       }
00120     }
00121   }
00122 }
00123 
00127 template <class TPSolver>
00128 void
00129 QCMDSystem<TPSolver>::fillMaps ( const vector< QCAtomIn > & v ) {
00130   int nbsbd = v.size() ; 
00131   if (nbsbd > 1) {
00132     for (int i = 0 ; i < nbsbd-1; ++i) {
00133       for (int j = i+1 ; j < nbsbd; ++j) {
00134         if (v[i].typeZone != QC_CORE && v[j].typeZone != QC_SHELL2) {
00135           domains[v[i].numDomain].addOverlpIndirectionInfo(v[j].numDomain, v[i].localNum, v[j].localNum,v[j].typeZone);
00136         }
00137         if (v[j].typeZone != QC_CORE && v[i].typeZone  != QC_SHELL2) {
00138           domains[v[j].numDomain].addOverlpIndirectionInfo(v[i].numDomain, v[j].localNum, v[i].localNum,v[i].typeZone);
00139           
00140         }
00141       }
00142     }
00143   }
00144 }
00145 
00149 template <class TPSolver>
00150 template <class TPSystem>
00151 void 
00152 QCMDSystem<TPSolver>::allocatePartitioner (QCPartitioner<TPSystem>*& partitioner, const QCGeneralData&  data){
00153   QC_TRACE_PART("BEGIN QCMDSystem<TPSolver>::allocatePartitioner ");
00157   QCPartitionerType pt = data.getPartitionerType();
00158   
00159   switch (pt) {
00160     
00161   case QC_LINEAR:
00162     partitioner = new QCLinearPartitioner<TPSystem>();
00163     break;
00164     
00165   case QC_KIMIKA:
00166     partitioner = new QCKimikaPartitioner<TPSystem>();
00167     break;
00168   default :
00169     std::cerr <<" Error on the name of the partitionner "<<pt<<". Cannot allocate the object. " << std::endl;
00170     exit(EXIT_FAILURE) ;
00171     
00172   }
00173   //
00174   QC_TRACE_PART("END QCMDSystem<TPSolver>::allocatePartitioner ");
00175 }
00176 
00177 template <class TPSolver>
00178 template <class TPSystem>
00179 void 
00180 QCMDSystem<TPSolver>::buildPartitioner (QCPartitioner<TPSystem>*& partitioner,
00181                                         const QCGeneralData&      data,
00182                                         const QCFiles&            files,
00183                                         const int *               charge,
00184                                         const int                 nbAtoms) {
00185   QC_TRACE_PART("BEGIN QCMDSystem<TPSolver>::buildPartitioner ");
00186   this->allocatePartitioner(partitioner,data) ;
00187   //
00188   partitioner->init(data,files,charge, nbAtoms);  // TO DEBUG A suuprimer
00189   //
00190   partitioner->setLoad(charge, nbAtoms);
00191   //
00192   QC_TRACE_PART("END QCMDSystem<TPSolver>::buildPartitioner ");
00193 
00194 }
00198 template <class TPSolver>
00199 template <class TPManager>
00200 void 
00201 QCMDSystem<TPSolver>::init (TPManager&    manager,  const string& path)
00202 {
00203   QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::init ");
00204 
00205   /***************************************************/
00206   QCPartitioner<TQCMDSystem> * partitioner = manager.getPartitioner();
00207 
00208   if (manager.getGeneralData().isReadPartitionFromFile()) { 
00209     if( partitioner ==NULL){
00210       this->allocatePartitioner(partitioner,const_cast<QCGeneralData&>(manager.getGeneralData()));
00211     }
00212     partitioner->readFromFile(*this, const_cast<QCGeneralData&>(manager.getGeneralData()),
00213                               manager.getFiles(),  path);
00214   } 
00215   else { 
00216     partitioner->initData(const_cast<QCGeneralData&>(manager.getGeneralData()), manager.getFiles());
00217     if(partitioner->getPartitionType() == QC_DIXON_PART  && 
00218        partitioner->getNbPartitions()!= this->getNbAA() ) {
00219        std::cerr << "Error For Dixon partitioning the number of partitions ("<<partitioner->getNbPartitions()<<")"
00220                  << " must be equal to the number of amino acids ("<<this->getNbAA() <<")"<<std::endl;
00221        exit(EXIT_FAILURE);
00222     }
00223     partitioner->partitioning(*this);
00224   }
00225   _partitionType = partitioner->getPartitionType() ;
00226   //
00228   //
00229   const QCGeneralData& data = manager.getGeneralData();
00230   QCIterator           workingSystem;
00231   if(data.getComputationType() != QC_TOOLS){
00232     for (workingSystem  = begin(); workingSystem != end();    workingSystem++) {
00233       if (data.getComputationType() == QC_ENERGY_DERIV) {
00234         workingSystem->allocateDisps();
00235       }
00236       workingSystem->sortOverlapments();
00237       workingSystem->setSystemNumbers(manager.getParameters());
00238       workingSystem->allocateMatrices(manager, 
00239                                       workingSystem->getNbAtomicOrbitals());
00240       workingSystem->fillWeight(manager.getGeneralData());
00241     }
00242   }
00243   QC_TRACE_INIT("END QCMDSystem<TPSolver>::init ");
00244   
00245 }
00249 template <class TPSolver>
00250 template <class TPManager>
00251 void 
00252 QCMDSystem<TPSolver>::completeHamiltonMatrices (TPManager& manager) {
00253 
00254   typedef typename TPManager::TModel         TModel;
00255   typedef typename TPManager::TModel::TParam TParam;
00256   //
00257   // On recupere avec des restricts pour le calcul les
00258   //  structures qui contiennent les donnees a manipuler
00259   //
00260   TModel&         QCRestrict model       = manager.getModel();
00261   const TParam *  QCRestrict params      = manager.getParameters();
00262   //
00263   //
00264   QCSymMatrix& QCRestrict hamiltonHAA = model.getSpWorkingAA();
00265   //
00266   // Les buffers du pipeline
00267   //
00268   const QCSubDomain ** contributers = new const QCSubDomain * [nbDomains];
00269   //
00270 #ifdef QC_VERBOSE
00271   std::cout << "* HAMILTON CONTRIBUTIONS *" << std::endl;
00272 #endif
00273   //  
00274   // En sequentiel les buffers sont de simples pointeurs
00275   //
00276   for (int i=0; i<nbDomains; ++i) {
00277     contributers[i] = &domains[i];
00278   }
00279   //    
00280   // On commence l'anneau de communication
00281   //
00282   for (int step=0; step < nbDomains-1; ++step) {
00283     
00284     // On simule la reception d'un domaine
00285     const QCSubDomain * recvd = contributers[nbDomains-1];
00286     
00287     // On avance dans le pipeline
00288     for (int i = nbDomains-1; i>0; --i) {
00289       contributers[i] = contributers[i-1];
00290     }
00291     contributers[0] = recvd;
00292     
00293 #ifdef QC_DEBUG_HAMILTON_CONTRIB_LEVEL2
00294     std::cout << std::endl << " *** STEP " << step << std::endl;
00295 #endif
00296     
00297     for (int i=0; i<nbDomains; ++i) {
00298       
00299 #ifdef QC_DEBUG_HAMILTON_CONTRIB_LEVEL2
00300       std::cout << "   * sd "  << domains[i].getId()  << " <= sd " << contributers[i]->getId() << std::endl;
00301 #endif
00302       
00303       domains[i].getHamiltonH().completeElems (// Les params
00304                                                params,
00305                                                // le domaine courant
00306                                                domains[i],
00307                                                // le domaine contribuant
00308                                                *contributers[i],
00309                                                // autres
00310                                                domains[i].getRepInterElecIntegrals(),
00311                                                hamiltonHAA, manager.getGeneralData(), step, nbDomains);
00312 
00313     }
00314   }
00315   
00316 #ifdef QC_OUTPUT_HAMILTON_CONTRIB 
00317   for (int i=0; i<nbDomains; ++i) {
00318     ostringstream osstr;
00319     osstr << "hamilton_mine_" << domains[i].getId()<<"_"; 
00320     domains[i].getHamiltonH().getMatrix().printInFile(osstr.str().c_str());
00321   }
00322 #endif
00323   delete [] contributers ;
00324 
00325 }
00326 
00330 template <class TPSolver>
00331 template <class TPManager>
00332 QCFloat 
00333 QCMDSystem<TPSolver>::completeFockMatricesAndElecEnergy (TPManager& manager, bool isFirstCall, int /* iter */) {
00334   //
00335   QC_TRACE("BEGIN QCMDSystem<TPSolver>::completeFockMatrices ");
00336   typedef typename TPManager::TModel         TModel;
00337   typedef typename TPManager::TModel::TParam TParam;
00338   
00343   TModel&         QCRestrict model  = manager.getModel();
00344   const TParam *  QCRestrict params = manager.getParameters();
00345 
00346   // Les parametrage de l'application
00347   const QCGeneralData& data = manager.getGeneralData();
00348 
00349   // matrices de Hamilton locales a cette methode
00350   QCSymMatrix& QCRestrict fockFAA    = model.getSpWorkingAA();
00351 
00352   // Les buffers du pipeline
00353   const QCSubDomain ** contributers  = new const QCSubDomain * [nbDomains];
00354   QCFloat           ** densityArrays = new QCFloat * [nbDomains];
00355   //
00356   // L'energie electronique sur tous les domaines
00357   QCFloat domainElecEnergy;
00358 
00359 
00360 #ifdef QC_VERBOSE
00361   std::cout << "* FOCK CONTRIBUTIONS *" << std::endl;
00362 #endif
00363   //
00364   // En memoire partagée les buffers sont de simples pointeurs
00365   for (int i=0; i < nbDomains; ++i) {
00366     contributers[i] = &domains[i];
00367     domains[i].initWeightedDensity(params);
00368     domains[i].fillDensityArray(params, densityArrays[i], true);
00369   }
00370   //  
00371   // On commence l'anneau de communication
00372   //
00373   for (int step=0; step < nbDomains-1; ++step) {
00374     
00375     // On simule la reception d'un domaine
00376     const QCSubDomain * recvdsbd = contributers[nbDomains-1];
00377     QCFloat           * recvdmat = densityArrays[nbDomains-1];
00378     
00379     // On avance dans le pipeline
00380     for (int i = nbDomains-1; i>0; --i) {
00381       contributers[i] = contributers[i-1];
00382       densityArrays[i] = densityArrays[i-1];
00383     }
00384 
00385     contributers[0]  = recvdsbd;
00386     densityArrays[0] = recvdmat;
00387 
00388 #ifdef QC_DEBUG_FOCK_CONTRIB_LEVEL2
00389     std::cout << std::endl << " *** STEP " << step << std::endl;
00390 #endif
00391 
00392     // Les contributions domaine par domaine
00393     for (int i=0; i < nbDomains; ++i) {
00394       
00395 #ifdef QC_DEBUG_FOCK_CONTRIB_LEVEL2
00396       std::cout << "   * sd " << domains[i].getId() << " <= sd " << contributers[i]->getId() 
00397            << std::endl;
00398 #endif
00399       
00400       domains[i].getFockF().completeElems (// Les params
00401                                            params,
00402                                            // le domaine courant
00403                                            domains[i],
00404                                            // le domaine contribuant
00405                                            *contributers[i], densityArrays[i],
00406                                            // autres
00407                                            domains[i].getRepInterElecIntegrals(),
00408                                            fockFAA, data, step, nbDomains, data.getIntgAcquisitionMethod(),
00409                                            isFirstCall);
00410       
00411     }
00412   }
00413 #if defined (QC_TRACE_IN_FILE)  
00414   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00415     qctrace[i].writeFockMatrix();
00416     }
00417 #endif
00418 
00419 #ifdef QC_OUTPUT_FOCK_CONTRIB
00420 
00421   for (int i=0; i < nbDomains; ++i) {
00422     ostringstream osstr;
00423     osstr << "fock_mine_" << domains[i].getId();
00424     domains[i].getFockF().getMatrix().printInFile(osstr.str().c_str());
00425   }
00426 
00427   std::cout << std::endl;
00428   
00429 #endif
00430   bool rewind =  (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE && !isFirstCall);
00431   //
00432   // On prepare le calcul de l'energie
00433   //
00434   QCEnergies::elecEnergy = QC_ZERO;
00435   //
00436   for (int i=0; i < nbDomains; ++i) {
00437     //
00438     // Liberation des buffers des blocs diagonaux de la matrice densite
00439     delete [] densityArrays[i];
00440     //
00441     // Retour au debut des fichiers d'intégrales de repultion
00442     if (rewind) {
00443       domains[i].getRepInterElecIntegrals().getIntgReader()->rewindFile();
00444     }
00445     //
00446     // Calcul de l'energie electronique sur chaque domaine
00447     QCIterator domainIter(this, i);
00448     domainElecEnergy = model.computeElecEnergy(domainIter, *domains[i].getWeightedDensityP());
00449 
00450 #ifdef VERBOSE_ENERGY
00451     std::cout << setprecision(20) << " * Electronic energy (sd " << i << ") = " << domainElecEnergy << std::endl;
00452 #endif
00453  #if defined (QC_TRACE_IN_FILE)  
00454     qctrace[i].writeVal("* Local ELECTRONIC ENERGY of my domain  = " ,domainElecEnergy);
00455 #endif
00456     QCEnergies::elecEnergy += domainElecEnergy;
00457   }
00458 
00459 #if defined (QC_TRACE_IN_FILE)  
00460   for (int i=0; i < nbDomains; ++i) {
00461     qctrace[i].writeVal("* Local ELECTRONIC ENERGY (ALL Domains in my MDDomain)  = " ,this->elecEnergy);
00462   }
00463 #endif
00464 
00465 #ifdef VERBOSE_ENERGY
00466   std::cout << setprecision(20)    << "* Total electronic energy  = " << QCEnergies::elecEnergy << std::endl;
00467 #endif
00468 #if defined (QC_TRACE_IN_FILE)  
00469   for (int i=0; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
00470     qctrace[i].writeVal("* Global ELECTRONIC ENERGY (all domains)  = " ,QCEnergies::elecEnergy);
00471   }
00472 #endif
00473   //
00474   // Free 
00475   //
00476   delete [] contributers;
00477   delete [] densityArrays;
00478 
00479   QC_TRACE("END   QCMDSystem<TPSolver>::completeFockMatrices ");
00480   return QCEnergies::elecEnergy;
00481 }
00482 
00486 template <class TPSolver>
00487 template <class TPMDSystem>
00488 void QCMDSystem<TPSolver>::adjustFermiEnergyNew (TPMDSystem& mdsystem, QCMemory& memory) {
00489   //
00490   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::adjustFermiEnergyNew ");
00491   // Les tableaux locaux
00492   QCFloat * orbitalOccupN;
00493   QCFloat evMin, evMax;
00494   QCFloat sumNE;
00495   QCWay way;
00496   int numSD;
00497   int * nbOM =   memory.template takeTmpMem<int>(nbDomains);
00498   QCTopEnergyLevel * topLevel =  memory.template takeTmpMem<QCTopEnergyLevel>(nbDomains);
00499   //
00500   // On met le vecteur d'occupation a 0
00501   //
00502   for (numSD = 0; numSD < nbDomains; ++numSD) {
00503     orbitalOccupN = domains[numSD].getOrbitalOccupN();
00504     nbOM[numSD]   = domains[numSD].getNbAtomicOrbitals();    
00505     memset(orbitalOccupN, 0, nbOM[numSD] * sizeof(QCFloat));
00506   }
00507   //
00508   // Adjustment
00509   //
00510   if (fermiFCall) {
00511     sumNE = QC_ZERO;
00512     way   = ASC;
00513     computeDichotomyBounds(mdsystem, nbOM, evMin, evMax);
00514     adjustFermiEnergyDichotomy(mdsystem, nbOM, topLevel, sumNE, evMin, evMax, way);
00515     } 
00516   else {
00517     //adjustFermiEnergyMonotone(mdsystem, nbOM, topLevel);
00518     adjustFermiEnergy(mdsystem, memory);
00519   }
00520   //
00521   fermiFCall = false;
00522   memory.template giveBackTmpMem<QCTopEnergyLevel>(topLevel);
00523   memory.template giveBackTmpMem<int>(nbOM);
00524   //
00525   // OPTIMAL DAMPING
00526   //
00527   QCFloat sumOccupEigenVal = 0.0;
00528   int nbOccupiedOM ;
00529   for (numSD = 0; numSD < this->nbDomains; ++numSD) {
00530     nbOccupiedOM                  = this->domains[numSD].getNbOccupiedOM();
00531     QCFloat * QCRestrict eigenVal = this->domains[numSD].getOrbitalOccupN();  
00532     //
00533     for (int i = 0 ; i < nbOccupiedOM ; ++i){
00534       sumOccupEigenVal += eigenVal[i] ;
00535     }    
00536   }
00537   //
00538   //
00539   QC_TRACE_ENER("END  QCMDSystem<TPSolver>::adjustFermiEnergyNew " << sumOccupEigenVal);
00540 }
00544 template <class TPSolver>
00545 template <class TPMDSystem>
00546 void QCMDSystem<TPSolver>::computeDichotomyBounds (TPMDSystem& system,
00547                                               int *nbOM, QCFloat& evMin, QCFloat& evMax) {
00548 
00549 
00550   int numSD, numOM;
00551   QCFloat * eigenVal;
00552   //
00553   evMin =  DBL_MAX;
00554   evMax = -DBL_MAX;
00555   
00556   for (numSD = 0; numSD < nbDomains; ++numSD) {
00557       
00558     eigenVal      = domains[numSD].getEigenVal();      
00559     nbOM[numSD]   = domains[numSD].getNbAtomicOrbitals();
00560     
00561     for (numOM = 0; numOM < nbOM[numSD]; ++numOM) {
00562       if (eigenVal[numOM] < evMin) {
00563         evMin = eigenVal[numOM];
00564       }
00565       if (eigenVal[numOM] > evMax) {
00566         evMax = eigenVal[numOM];
00567       }
00568     }
00569   }
00570   
00571   system.getGlobalMinMax(evMin, evMax);
00572 }
00573 
00577 template <class TPSolver>
00578 template <class TPMDSystem>
00579 void 
00580 QCMDSystem<TPSolver>::adjustFermiEnergyDichotomy(TPMDSystem&        system,
00581                                                  int *              nbOM,
00582                                                  QCTopEnergyLevel * topLevel,
00583                                                  QCFloat            sumNE,
00584                                                  QCFloat            evMin,
00585                                                  QCFloat            evMax,
00586                                                  QCWay              way) {
00587   
00588   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::adjustFermiEnergyDichotomy" );
00589   int loop, nbEV, numOM, numSD;
00590 
00591   QCFloat * bFactor;
00592   QCFloat * eigenVal;
00593   QCFloat * orbitalOccupN;
00594 
00595   //
00596   QCFloat dsumNE;
00597 
00598   // Dichotomy bounds
00599   QCFloat fMin, fMed, fMax;
00600   QCFloat fMedPrev;
00601   QCFloat gtEV, leEV;
00602   QCFloat buff[2];
00603   QCFloat diffNe;
00604 
00605 
00606   // Init bounds
00607   fMin = evMin;
00608   fMax = evMax;
00609   fMed = (fMin + fMax)*QC_HALF;
00610   
00611   // Fermi adjust loop
00612   const int MAX_LOOP = 500;
00613   loop = 0;
00614   while (loop < MAX_LOOP) {
00615 
00616     // Compute dsumNE
00617     dsumNE = QC_ZERO;
00618     nbEV   = 0;
00619     gtEV   = -DBL_MAX;
00620     leEV   =  DBL_MAX;
00621     
00622     bool first;
00623 
00624     for (numSD = 0; numSD < nbDomains; ++numSD) {
00625       
00626       eigenVal      = domains[numSD].getEigenVal();
00627       bFactor       = domains[numSD].getBFactor();
00628       nbOM[numSD]   = domains[numSD].getNbAtomicOrbitals();
00629       first         = true;
00630       
00631       for (numOM = 0; 
00632            (numOM < nbOM[numSD]) && (eigenVal[numOM]) <= fMax; 
00633            ++numOM) {
00634         
00635         if (eigenVal[numOM] < fMin) {
00636           continue;
00637         }
00638         if (eigenVal[numOM] < fMed && way == ASC) {
00639           dsumNE += bFactor[numOM];
00640         }
00641         if (eigenVal[numOM] >= fMed) {
00642           if (way == DESC) {
00643             dsumNE -= bFactor[numOM];
00644           }
00645           if (first) {
00646             topLevel[numSD].numOM = numOM;
00647             first = false;
00648           }
00649         }
00650         ++nbEV;
00651       }
00652       
00654       if (numOM > 0 && gtEV < eigenVal[numOM-1]) {
00655         gtEV = eigenVal[numOM-1];
00656         
00657       }
00658       if (numOM < nbOM[numSD] && leEV > eigenVal[numOM]) {
00659         leEV = eigenVal[numOM];
00660         
00661       }
00663     }
00664 
00665     buff[0] = dsumNE;
00666     buff[1] = nbEV;
00667     
00668     system.getTotalSumNbElecs(buff, 2);
00669 
00670     dsumNE = buff[0];
00671     nbEV   = static_cast<int>(buff[1]);
00672     
00673 
00674     sumNE   += QC_TWO*dsumNE;
00675     diffNe   = sumNE - nbElectrons;
00676     
00677 
00678 #ifdef QC_VERBOSE_FERMI_LEVEL2    
00679     std::cout << "** LOOP " << loop  << ": diffNE = " << diffNe << " , FEnergy = "
00680          << fMed << ", nbEV=" << nbEV << ", gtEV=" << gtEV << ", leEV=" << leEV  << std::endl;
00681 #endif
00682     //
00683     fMedPrev = fMed;
00684     if (diffNe > 0) {
00685       fMax = fMed;
00686       fMed = (fMed + fMin)*QC_HALF;
00687       way  = DESC;
00688 
00689     } else if (diffNe < 0) {
00690       fMin = fMed;
00691       fMed = (fMed + fMax)*QC_HALF;
00692       way  = ASC;
00693     }
00694     ++loop;
00695 
00696     if (nbEV == 1) {
00697       break;
00698     }
00699     
00700   }
00701 
00702   system.getGlobalMinMax(leEV, gtEV);
00703 
00704   // Set Fermi Energy
00705   if (way == ASC) {
00706     deltaFermiEnergy = FQCAbs(leEV - fermiEnergy);
00707     fermiEnergy = leEV;
00708 
00709   }
00710   if (way == DESC) {
00711     deltaFermiEnergy = FQCAbs(gtEV - fermiEnergy);
00712     fermiEnergy = gtEV;
00713   }
00714 
00715 #ifdef QC_VERBOSE_FERMI
00716   if (system.getRank() == 0) {
00717     std::cout << "PROC 0" << ": Fermi Energy = " << fermiEnergy << " (" << loop << " loops)" 
00718          << ", sumNE = " << sumNE << std::endl;
00719   }
00720 #endif
00721   //
00722   // Set orbitals occupation
00723   //
00724   for (numSD = 0; numSD < nbDomains; ++numSD) {
00725     eigenVal      = domains[numSD].getEigenVal();
00726     bFactor       = domains[numSD].getBFactor();
00727     orbitalOccupN = domains[numSD].getOrbitalOccupN();
00728     nbOM[numSD]   = domains[numSD].getNbAtomicOrbitals();
00729 
00730     for (numOM = 0;  numOM < nbOM[numSD] && eigenVal[numOM] < fermiEnergy; ++numOM) {
00731       orbitalOccupN[numOM] = QC_TWO;
00732     }
00733     domains[numSD].setNbFilledOccupiedOM(numOM, numOM);
00734     if (numOM < nbOM[numSD] && eigenVal[numOM] == fermiEnergy) {
00735       
00736       if (way == DESC) {
00737         orbitalOccupN[numOM] = QC_TWO + (nbElectrons - sumNE) / bFactor[numOM];
00738       
00739       } else if (way == ASC) {
00740         orbitalOccupN[numOM] = (nbElectrons - sumNE) / bFactor[numOM];
00741       }
00742         
00743       domains[numSD].setNbFilledOccupiedOM(numOM, numOM+1);
00744     
00745     } 
00746 
00747 
00748 #ifdef QC_VERBOSE_FERMI_LEVEL2
00749     std::cout << "sd" << domains[numSD].getId() << ": nbFilled = "
00750          << domains[numSD].getNbFilledOM() << ", nbOccupied = " << domains[numSD].getNbOccupiedOM()
00751          << std::endl;
00752 #endif
00753   }  
00754   
00755   QC_TRACE_ENER("END  QCMDSystem<TPSolver>::adjustFermiEnergyDichotomy" );
00756 }
00757 
00761 template <class TPSolver>
00762 template <class TPMDSystem>
00763 void 
00764 QCMDSystem<TPSolver>::adjustFermiEnergyMonotone(TPMDSystem&        system,
00765                                                 int *              nbOM,
00766                                                 QCTopEnergyLevel * topLevel) {
00767   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::adjustFermiEnergyMonotone   EMPTY"  );
00768   QC_TRACE_ENER("END QCMDSystem<TPSolver>::adjustFermiEnergyMonotone");
00769 }
00773 template <class TPSolver>
00774 template <class TPMDSystem>
00775 void 
00776 QCMDSystem<TPSolver>::adjustFermiEnergy (TPMDSystem& mdsystem, QCMemory& memory) {
00777   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::adjustFermiEnergy");
00778   QCFloat sumNbElecs;
00779   //
00780   // Les tableaux locaux
00781   QCFloat * eigenVal;
00782   QCFloat * bFactor;
00783   QCFloat * orbitalOccupN;
00784 
00785   // Les numeros
00786   int numOM, numSD;
00787   //
00788   QCTopEnergyLevel * topEnergyLevels = 
00789     memory.template takeTmpMem<QCTopEnergyLevel>(nbDomains);
00790   
00791   int * nbOM =  memory.template takeTmpMem<int>(nbDomains);
00792 
00793 
00794   sumNbElecs = QC_ZERO;
00795   for (numSD = 0; numSD < nbDomains; ++numSD) {
00796     numOM         = 0;
00797     eigenVal      = domains[numSD].getEigenVal();
00798     bFactor       = domains[numSD].getBFactor();
00799     orbitalOccupN = domains[numSD].getOrbitalOccupN();
00800     nbOM[numSD]   = domains[numSD].getNbAtomicOrbitals();
00801 
00802     // On met le vecteur d'occupation a 0
00803     memset(orbitalOccupN, 0, nbOM[numSD] * sizeof(QCFloat));
00804 
00805     while (numOM < nbOM[numSD] && eigenVal[numOM] < fermiEnergy) {
00806       sumNbElecs += bFactor[numOM];
00807 #ifdef QC_VERBOSE_FERMI_LEVEL2
00808 //      std::cout << "  -- numOM : "<< numOM << ", sumNbElecs: " << setprecision(6) << sumNbElecs << std::endl;
00809 #endif
00810       orbitalOccupN[numOM] = QC_TWO;
00811       ++numOM;
00812     }
00813     topEnergyLevels[numSD].setAttributes(numSD, numOM-1);
00814   }  
00815   //  
00816   // On multiplie par 2 le nb d'occupation de base des niveaux d'energie
00817   sumNbElecs *= QC_TWO;
00818   //
00819 #ifdef QC_VERBOSE_FERMI_LEVEL2
00820   std::cout << std::endl << "    -- Avant, sumNbElecs = " << sumNbElecs << std::endl;
00821   for (numSD = 0; numSD < nbDomains; ++numSD) {
00822     numOM = 0;
00823     orbitalOccupN = domains[numSD].getOrbitalOccupN();
00824     while (numOM < nbOM[numSD]) {
00825       std::cout << "    -- orbitalOccupN[" << numOM << "] = " << orbitalOccupN[numOM] << std::endl;
00826       ++numOM;
00827     }
00828   }  
00829 #endif
00830   mdsystem.getTotalSumNbElecs(&sumNbElecs, 1);
00831 
00832 
00833 #ifdef QC_VERBOSE_FERMI  
00834   //  if (mdsystem.getRank() == 0) {
00835   std::cout << std::endl  << "   --  sumNbElecs= " << sumNbElecs << "  nbElectrons= " << nbElectrons << std::endl;
00836     //  }
00837 #endif
00838   //
00839   if (sumNbElecs > nbElectrons) {
00840     mdsystem.removeLevels(sumNbElecs, topEnergyLevels);
00841   }
00842   if (sumNbElecs < nbElectrons) {
00843     mdsystem.addLevels(sumNbElecs, nbOM, topEnergyLevels);
00844   }
00845 
00846 #ifdef QC_VERBOSE_FERMI_LEVEL2
00847   std::cout << "   -- nbElecOfSystem = " << nbElectrons << "  fermiEnergy = " << fermiEnergy << std::endl;
00848   
00849   for (numSD = 0; numSD < nbDomains; ++numSD) {
00850     eigenVal      = domains[numSD].getEigenVal();
00851     bFactor       = domains[numSD].getBFactor();
00852     orbitalOccupN = domains[numSD].getOrbitalOccupN();
00853     
00854      numOM = 0;
00855      while (numOM < nbOM[numSD]) {
00856        std::cout << "  -- curEigenVal[" << numOM << "] = " <<  eigenVal[numOM] << " donc ";
00857        std::cout << "  -- curOrbitalOccupN[" << numOM << "] = " << orbitalOccupN[numOM] << " avec ";
00858        std::cout << "  -- curFactorB[" << numOM << "] = " << bFactor[numOM] << std::endl;      
00859        ++numOM;
00860      }
00861     std::cout << "  -- sd" << numSD << ": nbFilled = "   << domains[numSD].getNbFilledOM()
00862          << ", nbOccupied = " << domains[numSD].getNbOccupiedOM() << std::endl;
00863   }
00864 #endif
00865   //
00866   memory.template giveBackTmpMem<int>(nbOM);
00867   memory.template giveBackTmpMem<QCTopEnergyLevel>(topEnergyLevels);
00868   QC_TRACE_ENER("END   QCMDSystem<TPSolver>::adjustFermiEnergy");
00869 }
00873 template <class TPSolver>
00874 template <class TPManager>
00875 void QCMDSystem<TPSolver>::completeDensityMatrices (TPManager& manager) {
00876   QC_TRACE_DENSITY("BEGIN QCMDSystem<TPSolver>::completeDensityMatrices");
00877 
00878   typedef typename TPManager::TModel         TModel;
00879   typedef typename TPManager::TModel::TParam TParam;
00880   
00885   TModel&         QCRestrict model  = manager.getModel();
00886   const TParam *  QCRestrict params = manager.getParameters();
00887   
00888 
00889   // Les matrices densite
00890   QCDensityGtr<QCSymMatrix> * QCRestrict densityP;
00891   QCSymMatrix * QCRestrict weightedDensity;
00892   QCSymMatrix * QCRestrict remWeightedDensity;
00893 
00894   // matrices de Hamilton locales a cette methode
00895   QCMatrix&    QCRestrict wDensityAB = model.getSpWorkingAB();
00896   QCMatrix&    QCRestrict interPAB   = model.getSpWorkingAB2();
00897   QCSymMatrix& QCRestrict wDensityAA = model.getSpWorkingAA();
00898 
00899 
00900 #ifdef QC_VERBOSE_DENSITY
00901   std::cout << "*       DENSITY CONTRIBUTIONS *" << std::endl;
00902 #endif
00903   
00904   
00905   for (int i = 0; i < nbDomains; ++i) {
00906 
00907     densityP        = &domains[i].getDensityP();
00908     weightedDensity = &domains[i].getWeightedDensityP()->getMatrix();
00909     
00910     // On copie la matrice densite ponderee dans
00911     // la matrice densite
00912     densityP->getMatrix().copy(*weightedDensity);
00913 
00914 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
00915     std::cout << "   BEGIN Contribution for domain "   << domains[i].getId()   << " (local number " << i    << std::endl;
00916 #endif
00917     for (int j = 0; j < domains[i].getNbOverlpSD(); ++j) {
00918 
00919       int rId = domains[i].getRemDomainId(j);
00920       remWeightedDensity = &domains[rId].getWeightedDensityP()->getMatrix();
00921 
00922 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
00923       std::cout << "               Contribution from  sd " << domains[rId].getId()    << std::endl;
00924 #endif
00925       
00926       densityP->mergeDensity(params,domains[i],domains[rId],*remWeightedDensity,j,wDensityAB,interPAB, wDensityAA);
00927     }
00928 #ifdef QC_DEBUG_DENSITY_CONTRIB_LEVEL2
00929     std::cout << "   END Contribution for domain "   << domains[i].getId()   << " (local number " << i    << std::endl;
00930 #endif
00931 
00932 #ifdef QC_OUTPUT_DENSITY_CONTRIB
00933     ostringstream osstr;
00934     osstr << "merged_density_mine_" << domains[i].getId();
00935     densityP->getMatrix().printInFile(osstr.str().c_str());
00936 #endif
00937   }
00938   QC_TRACE_DENSITY("END   QCMDSystem<TPSolver>::completeDensityMatrices");
00939 
00940 }
00944 template <class TPSolver>
00945 template <class TPParam>
00946 QCFloat QCMDSystem<TPSolver>::applyOptimalDamping (const TPParam * QCRestrict params) {
00947   //
00948   QC_TRACE_SCF("BEGIN QCMDSystem<TPSolver>::applyOptimalDamping");
00949   //Variables pour l'Optimal damping.
00950   QCFloat globalTraceFP1_PP1OfSD, globalTraceF_PP1OfSD;
00951   QCFloat globalTraceFP1_POfSD, globalTraceF_POfSD;
00952   //
00953   QCFloat lambdaOptDamp;
00954   //
00955   globalTraceFP1_PP1OfSD = QC_ZERO;
00956   globalTraceF_PP1OfSD   = QC_ZERO;
00957   globalTraceFP1_POfSD   = QC_ZERO;
00958   globalTraceF_POfSD     = QC_ZERO;
00959   //
00960   for (int i=0; i < nbDomains; ++i) {
00961     globalTraceFP1_PP1OfSD += 
00962       domains[i].traceProductCoreShell1(params, domains[i].getFockF().getMatrix(),
00963                                         domains[i].getDensityP().getMatrix());
00964     
00965     globalTraceF_PP1OfSD   += 
00966       domains[i].traceProductCoreShell1(params, domains[i].getFockFTild()->getMatrix(),
00967                                         domains[i].getDensityP().getMatrix());
00968     
00969     globalTraceFP1_POfSD   += 
00970       domains[i].traceProductCoreShell1(params, domains[i].getFockF().getMatrix(),
00971                                         domains[i].getDensityPTild()->getMatrix());
00972     
00973     globalTraceF_POfSD     += 
00974       domains[i].traceProductCoreShell1(params, domains[i].getFockFTild()->getMatrix(),
00975                                         domains[i].getDensityPTild()->getMatrix());
00976   }
00977   
00978   lambdaOptDamp = calculateLambdaOptDamp(globalTraceFP1_PP1OfSD, globalTraceF_PP1OfSD,
00979                                          globalTraceFP1_POfSD, globalTraceF_POfSD);
00980   
00981   for (int i=0; i < nbDomains; ++i) {
00982     domains[i].modifyMatricesAndEnergy(lambdaOptDamp);
00983   }
00984   QC_TRACE_SCF("END QCMDSystem<TPSolver>::applyOptimalDamping");
00985   return lambdaOptDamp ;
00986 }
00987 
00991 template <class TPSolver>
00992 QCFloat QCMDSystem<TPSolver>::calculateLambdaOptDamp (const QCFloat& globalTraceFP1_PP1OfSD,
00993                                                       const QCFloat& globalTraceF_PP1OfSD,
00994                                                       const QCFloat& globalTraceFP1_POfSD,
00995                                                       const QCFloat& globalTraceF_POfSD) {
00996   //
00997   QC_TRACE_SCF("BEGIN QCMDSystem<TPSolver>::calculateLambdaOptDamp");
00998   QCFloat sOptDamp, cOptDamp, lambdaOptDamp;
00999   //
01000   // On n'a pas de facteur 2 ici car on est 
01001   //  en semi-empirique ou P = 2P (modele quantique classique)
01002   //
01003   sOptDamp = globalTraceF_PP1OfSD - globalTraceF_POfSD;
01004   //
01005   // On a un facteur 1/2 ici car on est en 
01006   // semi-empirique ou P = 2P(modele quantique classique)
01007   //
01008   cOptDamp      = QC_HALF * (globalTraceFP1_PP1OfSD - globalTraceFP1_POfSD - sOptDamp);
01009   lambdaOptDamp = QC_ONE;
01010   //
01011   // Attention, on protege la division par zero.
01012   // lambdaOptDamp = min { -sOptDamp / QCTwo * cOptDamp), 1}
01013   //
01014   if (cOptDamp > QC_ZERO && -sOptDamp < (QC_TWO * cOptDamp) ) {
01015     lambdaOptDamp = -sOptDamp / (QC_TWO * cOptDamp);
01016   }
01017   QC_TRACE_SCF("END QCMDSystem<TPSolver>::calculateLambdaOptDamp : " << lambdaOptDamp);
01018   //
01019   return lambdaOptDamp; //Enfin, on peut renvoyer lambdaOptDamp.
01020 }
01021 
01025 template <class TPSolver>
01026 bool QCMDSystem<TPSolver>::testConvergence (QCFloat & error, QCFloat threshold) {
01027   
01028   bool    convReached = true;
01029   QCFloat errorL;
01030   //
01031   errorL  = error  = QC_ZERO;
01032   for (int i = 0; i < nbDomains; ++i) {
01033     error   = std::max(error,
01034                         domains[i].getDensityP().deltaElemMax(domains[i].getDensityPnm1()));
01035     convReached = convReached && (error <= threshold);
01036 //     std::cout << "    Dom " << i 
01037 //            << " D error "<< domains[i].getDensityP().deltaElemMax(domains[i].getDensityPnm1()) <<std::endl;
01038   }
01039 
01040   return convReached;
01041 }
01042 
01046 template <class TPSolver>
01047 QCFloat
01048 QCMDSystem<TPSolver>::getTotalSumTrace (void) {
01049   QCFloat sumTrace = 0;
01050 
01051   for (int i=0; i < nbDomains; ++i) {
01052     sumTrace += domains[i].getWeightedDensityP()->trace();
01053   }
01054 
01055   return sumTrace;
01056 }
01060 template <class TPSolver>
01061 void
01062 QCMDSystem<TPSolver>::removeLevels (QCFloat& sumNbElecs,
01063                                     QCTopEnergyLevel * topEnergyLevels) {
01064   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::removeLevels -- Adjust Fermi Energy ");
01065   QC_TRACE_ENER("         nbDomains " << nbDomains  );
01066 
01067   //
01068   // On change toujours les elements de la derniere case
01069   //         car tri par ordre croissant
01070   //
01071   int& numSDToChange = topEnergyLevels[nbDomains-1].numSD;
01072   int& numOMToChange = topEnergyLevels[nbDomains-1].numOM;
01073   QCFloat& energyLevelToChange = topEnergyLevels[nbDomains-1].energyLevel;
01074   int numSD, numOfHomoOfSD;
01075 
01076   bool rightNumberOfElec = false;
01077 
01078   // On met dans topEnergyLevels les energies des homos.
01079   for (numSD = 0; numSD < nbDomains; ++numSD) {
01080     numOfHomoOfSD = topEnergyLevels[numSD].numOM;
01081     if (numOfHomoOfSD >= 0) {
01082       topEnergyLevels[numSD].energyLevel = domains[numSD].getEigenVal(numOfHomoOfSD);
01083     } 
01084     else {
01085       // On ne peut plus enlever d'electrons sur ce sous-domaine.
01086       topEnergyLevels[numSD].energyLevel = -DBL_MAX;
01087     }
01088   }
01089   //
01090   int loop = 0, choice = -1;
01091   while (!rightNumberOfElec) {
01092     
01093     // On trie les niveaux par ordre croissant.
01094     stable_sort(topEnergyLevels, topEnergyLevels + nbDomains);
01095 
01096     // On enleve le nb de du niveau le + haut au niveau au total.
01097     if ( (sumNbElecs - QC_TWO * domains[numSDToChange].getBFactor(numOMToChange)) >= nbElectrons) {
01098       choice = 0 ;
01099       sumNbElecs -= QC_TWO * domains[numSDToChange].getBFactor(numOMToChange);
01100       domains[numSDToChange].setOrbitalOccupN(numOMToChange, QC_ZERO);
01101 
01102       // On le remplace par le niveau le plus haut du meme sous domaine.
01103       --numOMToChange;
01104       if (numOMToChange >= 0) {
01105         choice = 1 ;
01106         energyLevelToChange = domains[numSDToChange].getEigenVal(numOMToChange);
01107       
01108       } else {
01109         // On ne peut plus enlever d electrons sur ce sous-domaine.
01110         choice = 2 ;
01111         energyLevelToChange = -DBL_MAX;
01112       }
01113     
01114     } 
01115     else {
01116       // On diminue le nombre d'occupation de l'homo.
01117       choice =3 ;
01118       domains[numSDToChange].setOrbitalOccupN(numOMToChange,  QC_TWO +
01119                                               (nbElectrons - sumNbElecs) / domains[numSDToChange].getBFactor(numOMToChange));
01120       rightNumberOfElec = true;
01121     }
01122     
01123     QC_TRACE_ENER("     loop : " << loop <<"        choice " << choice  );
01124     ++loop;
01125   }
01126   //
01127   // On repositionne l'energie de Fermi. Ici l'energie de Fermi est l'energie 
01128   // de la Homo
01129   fermiEnergy = energyLevelToChange;
01130   //
01131   QC_TRACE_ENER("        (" << loop << " loops) fermi energy = " << fermiEnergy );
01132   //
01133   // On affecte nbFilledOM et nbOccupiedOM du domaine
01134   //
01135   // Dans topEnergyLevels, on a les homos.
01136   for (numSD = 0; numSD < nbDomains; ++numSD) {
01137 
01138     // Le nb d'orbitale occupe est egal au nombre d'orbitales remplies.
01139     
01140     // Le sous-domaine est vide.
01141     if (topEnergyLevels[numSD].energyLevel == -DBL_MAX) {
01142       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(0, 0);
01143     }
01144     
01145     // L'homo du SD est remplie.
01146     else if (domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(topEnergyLevels[numSD].numOM)
01147              == QC_TWO) {
01148 
01149       // La derniere orbitale est completement remplie.
01150       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM + 1,
01151                               topEnergyLevels[numSD].numOM + 1);
01152     }
01153     
01154     // L'homo du SD est occupee.
01155     else {
01156 
01157       // La derniere orbitale n'est pas completement remplie.
01158       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM,
01159                               topEnergyLevels[numSD].numOM + 1);
01160     }
01161   }
01162   QC_TRACE_ENER("END   QCMDSystem<TPSolver>::removeLevels -- Adjust Fermi Energy");
01163 }
01164 
01168 template <class TPSolver>
01169 void QCMDSystem<TPSolver>::addLevels (QCFloat& sumNbElecs, const int * nbOM,
01170                                  QCTopEnergyLevel * topEnergyLevels) {
01171 
01172   QC_TRACE_ENER("BEGIN QCMDSystem<TPSolver>::addLevels" );  
01173   //
01174   // On change toujours les elements de la premiere case
01175   // car tri par ordre croissant
01176   int&     numSDToChange       = topEnergyLevels[0].numSD; 
01177   int&     numOMToChange       = topEnergyLevels[0].numOM;
01178   QCFloat& energyLevelToChange = topEnergyLevels[0].energyLevel;
01179   
01180   int numSD;
01181   int numOfLumoOfSD;
01182 
01183   bool rightNumberOfElec = false;
01184 
01185   // On met dans topEnergyLevels les energies des lumos.
01186   for (numSD = 0; numSD < nbDomains; ++numSD) {
01187     numOfLumoOfSD = topEnergyLevels[numSD].numOM + 1;
01188     if (numOfLumoOfSD < nbOM[numSD]) {
01189       topEnergyLevels[numSD].numOM = numOfLumoOfSD;
01190       topEnergyLevels[numSD].energyLevel = domains[numSD].getEigenVal(numOfLumoOfSD);
01191       
01192     } else {
01193       // On ne peut plus enlever d'electrons sur ce sous-domaine.
01194       topEnergyLevels[numSD].energyLevel = DBL_MAX;
01195     }
01196   }
01197   //
01198   int loop = 0;
01199   while (!rightNumberOfElec) {
01200     
01201     // On trie les niveaux par ordre croissant.
01202     stable_sort(topEnergyLevels, topEnergyLevels + nbDomains);
01203   
01204     // On ajoute le nb d'e du niveau le + haut au niveau au total.
01205     if ( (sumNbElecs + QC_TWO * domains[numSDToChange].getBFactor(numOMToChange)) <=
01206          nbElectrons) {
01207 
01208       sumNbElecs += QC_TWO * domains[numSDToChange].getBFactor(numOMToChange);
01209       domains[numSDToChange].setOrbitalOccupN(numOMToChange, QC_TWO);
01210 
01211       // On le remplace par le niveau le plus haut du meme sous domaine.
01212       numOMToChange++;
01213       if (numOMToChange < nbOM[numSDToChange]) {
01214         energyLevelToChange = domains[numSDToChange].getEigenVal(numOMToChange);
01215       
01216       } else {
01217         // On ne peut plus enlever d electrons sur ce sous-domaine.
01218         energyLevelToChange = DBL_MAX;
01219       }
01220     
01221     } else {
01222       // On augmente le nombre d'occupation de la lumo.
01223       domains[numSDToChange].setOrbitalOccupN(numOMToChange, 
01224                                               (nbElectrons - sumNbElecs) / domains[numSDToChange].getBFactor(numOMToChange));
01225      
01226       // La lumo est la suivante
01227       numOMToChange++;
01228       rightNumberOfElec = true;
01229     }
01230 
01231     ++loop;
01232   }
01233   
01234   // On repositionne l'energie de Fermi. Ici l'energie de Fermi est l'energie 
01235   // de la Lumo
01236   fermiEnergy = energyLevelToChange;
01237   //
01238   QC_TRACE_ENER("       (" << loop  << " loops) fermi energy = " << fermiEnergy );  
01239   //
01240   // On affecte nbFilledOM et nbOccupiedOM du domaine
01241   //    Dans topEnergyLevels, on a les lumos.
01242   //
01243   for (numSD = 0; numSD < nbDomains; ++numSD) {
01244     //
01245     // Le nb d'orbitale occupe est egal au nombre d'orbitales remplies.
01246     // Le sous-domaine est rempli
01247     if (topEnergyLevels[numSD].energyLevel == DBL_MAX) {
01248       domains[topEnergyLevels[numSD].numSD].
01249         setNbFilledOccupiedOM(nbOM[topEnergyLevels[numSD].numSD], 
01250                               nbOM[topEnergyLevels[numSD].numSD]);
01251     }
01252 
01253     // Le sous-domaine est vide
01254     else if (domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(0) == QC_ZERO) {
01255       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(0, 0);
01256     }
01257     // La lumo du SD est remplie
01258     else if (domains[topEnergyLevels[numSD].numSD].getOrbitalOccupN(topEnergyLevels[numSD].numOM-1)
01259              == QC_TWO) {
01260 
01261       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM,topEnergyLevels[numSD].numOM);
01262     }
01263     
01264     // La lumo du SD est juste occupee
01265     else {
01266   
01267       // La derniere orbitale n'est pas completement remplie.
01268       domains[topEnergyLevels[numSD].numSD].setNbFilledOccupiedOM(topEnergyLevels[numSD].numOM - 1,topEnergyLevels[numSD].numOM);
01269     } 
01270   } 
01271   QC_TRACE_ENER("END   QCMDSystem<TPSolver>::addLevels" );    
01272 }
01273 
01277 template <class TPSolver>
01278 ostream& operator << (ostream& out, const QCMDSystem<TPSolver>& system) {
01279   //
01280   out << std::endl << ">>>>> ROOT SYSTEM <<<<<"   << std::endl
01281       << std::endl << system.getRootSystem() << std::endl;
01282   for (int i=0; i < system.getNbDomains(); ++i) {
01283     out << std::endl << ">>>>> DOMAIN " << system.getSubDomain(i).getId() 
01284         << " <<<<<"   << std::endl << std::endl << system.getSubDomain(i);
01285   }
01286   return out;
01287 }
01288 
01289 template <class TPSolver>
01290 QCFloat  
01291 QCMDSystem<TPSolver>::computeSumOfEigenValOccupied (){
01292   QCFloat t = 0.0, tloc=0.0;
01293   QCFloat * eigenVal = NULL;
01294   for (int i = 0 ; i < QCMDSystem<TPSolver>::nbDomains; ++i) {
01295         eigenVal      = domains[i].getEigenVal();      
01296         for (int j = 0 ; j < domains[i].getNbOccupiedOM(); ++j) {
01297           tloc += eigenVal[j] ;
01298         }
01299   }
01300   return t ;
01301 }
01302 
01305 // template <class TPSolver>
01306 // template <class TPManager>
01307 // void QCMDSystem<TPSolver>::void computeElecEnergy(TPManager& manager){
01308 //   //
01309 //   typedef typename TPManager::TModel         TModel;
01310 //   TModel&         QCRestrict model  = manager.getModel();
01311 //   QCIterator           workingSystem;
01312 
01313 //   //
01314 //   // Boucle sur les SD
01315 //   //
01316 //   for (workingSystem  = this->begin(); workingSystem != this->end();    workingSystem++) {
01317 //     model.computeElecEnergy(current,*(current->getWeightedDensityP()));
01318 //   }
01319 // }
01320 template <class TPSolver>
01321 void 
01322 QCMDSystem<TPSolver>::readDensityFromFileAscii(const QCFiles & files,  const string& path){
01323   QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::readDensityOnFileAscii   "<<files.getDensityFile());
01324   //  
01325   std::ifstream data ;
01326   std::string fileName , type ;
01327   //  fileName = path + "/" + files.getDensityFile() ;
01328   fileName = files.getDensityFile() ;
01329   data.open(fileName.c_str());
01330   //
01331   if (!data){
01332     std::cerr << "Error to open density Matrix in File " << fileName <<std::endl; 
01333     exit(EXIT_FAILURE) ;
01334   }
01335   //
01336   int numberOfSD ;
01337   data >> type ; 
01338   if (type ==  "BASIC-SHARED"){
01339     data >> numberOfSD ;
01340     if(this->nbDomains != numberOfSD){
01341       std::cerr << "The number of subdomain is different between the current partitionning and this in density Matrix in File."
01342                 <<std::endl
01343                 << " Must be the same " <<std::endl;
01344       exit(EXIT_FAILURE) ;
01345     }   
01346     this->readSharedDensityFromFileAscii(data,fileName) ;
01347   }
01348   else if (type ==  "BASIC-PARA"){
01349     int numberOfProc ;
01350     data >> numberOfProc >> numberOfSD ;
01351      if(this->nbDomains != numberOfSD){
01352       std::cerr << "The number of subdomain is different between the current partitionning and this in density Matrix in File."
01353                 <<std::endl
01354                 << " Must be the same " <<std::endl;
01355       exit(EXIT_FAILURE) ;
01356     }   
01357    std::cerr << "NOT YET IMPLEMENTED for BASIC-PARA " <<std::endl;
01358    exit(EXIT_FAILURE) ;
01359    this->readDistDensityFromFileAscii(data,numberOfProc) ;
01360   }
01361   else {
01362     std::cerr << "Bad type ("<<type<<") of storage in density Matrix in File " << fileName <<std::endl; 
01363     exit(EXIT_FAILURE) ;
01364   }
01365   QC_TRACE_INIT("END   QCMDSystem<TPSolver>::readDensityOnFileAscii   ");
01366 }
01367 template <class TPSolver>
01368 void 
01369 QCMDSystem<TPSolver>::readSharedDensityFromFileAscii(std::ifstream& data, const std::string &fileName){
01370   QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::readSharedDensityOnFileAscii   ");
01371   //
01372   std::string type ;
01373   int iloc, jloc, iglob, jglob, flag, nsd, nbElements,dim, bidon;
01374   QCFloat val;
01375   for (int k = 0 ; k < this->nbDomains; ++k) {
01376     data >> type  >>  bidon >> nsd >> dim >> nbElements ;
01377     QCSymMatrix& QCRestrict  D = this->domains[nsd].getDensityP().getMatrix();
01378     if ( dim != D.getDim() ){
01379       std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
01380                 << " Dimension in File : " << dim <<"  Dimension given by the model" << D.getDim() <<std::endl ;
01381       exit (EXIT_FAILURE) ;
01382     }
01383     for (int e = 0; e < nbElements ; ++e ){
01384       data >> iloc >> jloc >> iglob >> jglob >> flag  >>val ;
01385       if(flag == 0){ continue; }
01386       if (iloc < jloc ){ 
01387         std::cerr <<" Error we must read a Symetric Matrix " <<std::endl ;
01388         exit(EXIT_FAILURE) ;
01389       }
01390       D[iloc][jloc] = val ;
01391     }
01392   }
01393   QC_TRACE_INIT("END   QCMDSystem<TPSolver>::readSharedDensityOnFileAscii   ");
01394 }
01395 
01396 template <class TPSolver>
01397 void 
01398 QCMDSystem<TPSolver>::readDistDensityFromFileAscii(std::ifstream& data, const int nProc){
01399   QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::readDistDensityOnFileAscii   ");
01400   //
01401   std::ifstream dataD;
01402   std::string fileName , type ;
01403   int iloc, jloc, iglob, jglob, flag, nsd, nsdloc, nbElements,dim,numberOfSD;
01404   QCFloat val;
01405   for(int p = 0 ; p < nProc; ++p){
01406     data >> fileName ;
01407     dataD.open(fileName.c_str());
01408     if(!dataD){
01409       std::cerr << "Error to open density Matrix in File " << fileName <<std::endl; 
01410       exit(EXIT_FAILURE) ;    }
01411     dataD >> type >> numberOfSD;
01412     //
01413     for (int k = 0 ; k < numberOfSD; ++k) {
01414       dataD >> type  >>  nsdloc >>nsd >> dim >> nbElements ;
01415       QCSymMatrix& QCRestrict  D = this->domains[nsd].getDensityP().getMatrix();
01416       if ( dim != D.getDim() ){
01417         std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
01418                   << " Dimension in File : " << dim <<"  Dimension given by the model" << D.getDim() <<std::endl ;
01419         exit (EXIT_FAILURE) ; }
01420       for (int e = 0; e < nbElements ; ++e ){
01421         data >> iloc >> jloc >> iglob >> jglob >> flag  >>val ;
01422         if(flag == 0){ continue; }
01423         if (iloc < jloc ){ 
01424           std::cerr <<" Error we must read a Symetric Matrix " <<std::endl ;
01425           exit(EXIT_FAILURE) ;
01426         }
01427         D[iloc][jloc] = val ;
01428       }
01429     }
01430   }
01431   QC_TRACE_INIT("END   QCMDSystem<TPSolver>::readDistDensityOnFileAscii   ");
01432 }
01433 template <class TPSolver>
01434 template <class TPManager>
01435 void 
01436 QCMDSystem<TPSolver>::writeDensityOnFileAscii(TPManager& manager,  const std::string& path){
01437 // template <class TPSolver>
01438 // void 
01439 // QCMDSystem<TPSolver>::writeDensityOnFileAscii(const QCFiles & files,  const string& path){
01440   QC_TRACE_OUT("BEGIN QCMDSystem<TPSolver>::writeDensityOnFileAscii   ");
01441   //
01442   std::ofstream out;
01443   std::string fileName;
01444   //
01445   fileName = path + "/" + manager.getFiles().getResultFile() + "-density.ascii";
01446   //
01447   out.open(fileName.c_str());
01448   //
01449   out.setf(std::ios_base::scientific) ;   out.precision(20) ;
01450   //
01451   // ici il faut faire une boucle sur les SD
01452   //
01453   int flag = 0 ;
01454   //
01455   //  Il faut construire la correspondance entre l'orbitale locale et l'orbitale globale
01456   // Les parametres
01457   const typename TPManager::TParam *QCRestrict parameters =   manager.getParameters();
01458  // Le parametre courant
01459   const typename TPManager::TParam *QCRestrict parameterA;
01460   const typename TPManager::TParam *QCRestrict parameterB;
01461   //
01462   int firstAOofA, firstAOofGA, firstAOofSucA, firstAOofB, firstAOofGB, firstAOofSucB;
01463   //
01464   out  << "BASIC-SHARED   "<< this->nbDomains << std::endl;
01465   for (int k = 0 ; k < this->nbDomains; ++k) { // loop on subdomain
01466     //
01467     //std::cout << this->domains[k] <<std::endl;
01468     //
01469     //
01470     QCSymMatrix& QCRestrict  D = this->domains[k].getDensityP().getMatrix();
01471     //
01472     //    il faut construire le flag
01473     //
01474     QCSubDomainZone  * zone         = this->domains[k].getZone() ;
01475     const int        * globalIndex  = this->domains[k].getGlobalIdxs() ;
01476     QCSubDomainZone  zonei, zonej ;
01477     QCPartitionType  partType       = manager.getPartitioner()->getPartitionType() ;
01478     int cptA, cptB ;
01479     //
01480     // Compute the number of elements to save
01481     //
01482     int numElts0 = 0 ;
01483     for (int i = 0; i <  this->domains[k].getNbAtoms(); ++i) {
01484     
01485       parameterA    = &parameters[ this->domains[k].getParamIndexAt(i)];
01486       firstAOofA    =  this->domains[k].getFirstAOAt(i);
01487       firstAOofSucA = firstAOofA + parameterA->getNbAO();
01488       //
01489       zonei = zone[i] ;
01490       for (int j = 0; j <= i; ++j) {
01491         parameterB    = &parameters[ this->domains[k].getParamIndexAt(j)]; 
01492         firstAOofB    =  this->domains[k].getFirstAOAt(j); 
01493         firstAOofSucB = firstAOofB + parameterB->getNbAO();
01494         //
01495         zonej = zone[j] ;
01496         flag = computeFlag(zonei,zonej,partType); 
01497         //
01498         if( flag == 0) {continue ;}
01499         for (int mu = firstAOofA; mu < firstAOofSucA; ++mu) {
01500           for (int nu = firstAOofB; nu < firstAOofSucB; ++nu) {
01501             if( mu < nu) { break;}
01502             ++numElts0 ;
01503           }
01504         } 
01505       }
01506     }
01507     out  << "BASIC  "<< k <<"   "<< k <<"   "<<  this->domains[k].getGlobalIdx(k) <<"   "
01508          << D.getDim() << "  "<< numElts0 << std::endl; 
01509     for (int i = 0; i <  this->domains[k].getNbAtoms(); ++i) {
01510     
01511       parameterA    = &parameters[ this->domains[k].getParamIndexAt(i)];
01512       firstAOofA    =  this->domains[k].getFirstAOAt(i);
01513       firstAOofGA   =  QCSystem::getFirstAOAt(globalIndex[i]); // accede au systèmeglobal
01514       firstAOofSucA = firstAOofA + parameterA->getNbAO();
01515       //
01516       zonei = zone[i] ;
01517       for (int j = 0; j <= i; ++j) {
01518         parameterB    = &parameters[ this->domains[k].getParamIndexAt(j)]; 
01519         firstAOofB    =  this->domains[k].getFirstAOAt(j); 
01520         firstAOofGB   =  QCSystem::getFirstAOAt(globalIndex[j]);
01521         firstAOofSucB = firstAOofB + parameterB->getNbAO();
01522         //
01523         zonej = zone[j] ;
01524         flag = computeFlag(zonei,zonej,partType); 
01525         //
01526         if( flag == 0) {continue ;}
01527         cptA = firstAOofGA - 1 ;
01528         for (int mu = firstAOofA; mu < firstAOofSucA; ++mu) {
01529           ++cptA ;  cptB = firstAOofGB - 1 ;
01530           for (int nu = firstAOofB; nu < firstAOofSucB; ++nu) {
01531             ++cptB ;
01532             if( mu < nu) { break;;}
01533             if(cptA < cptB){ //The global density matrix is also symtric
01534               out << "  "<< mu << "  "<< nu << " "<< cptB << "  " << cptA <<  "  " << flag 
01535                   << "  " << D[mu][nu] << std::endl;
01536             }
01537             else{
01538               out << "  "<< mu << "  "<< nu << " "<< cptA << "  " << cptB <<  "  " << flag << "  " 
01539                   << D[mu][nu] << std::endl;
01540             }
01541           }
01542         } 
01543       }
01544     }
01545   }
01546   // 
01547   out.close();
01548   QC_TRACE_OUT("END   QCMDSystem<TPSolver>::writeDensityOnFileAscii   " ) ;
01549 }
01550 template <class TPSolver>
01551 template <class TPManager>
01552 void 
01553 QCMDSystem<TPSolver>::writeLocalDensityOnFilesAscii(TPManager& manager,  const std::string& path){
01554   QC_TRACE_OUT("BEGIN QCMDSystem<TPSolver>::writeLocalDensityOnFilesAscii   ");
01555   //
01556   std::ofstream out;
01557   std::string fileName;
01558   //
01559   //  Il faut construire la correspondance entre l'orbitale locale et l'orbitale globale
01560   // Les parametres
01561   const typename TPManager::TParam *QCRestrict parameters =   manager.getParameters();
01562  // Le parametre courant
01563   const typename TPManager::TParam *QCRestrict parameterA;
01564   const typename TPManager::TParam *QCRestrict parameterB;
01565   //
01566   int firstAOofA, firstAOofGA, firstAOofSucA, firstAOofB, firstAOofGB, firstAOofSucB;
01567   //
01568   out  << "BASIC-SHARED   "<< this->nbDomains << std::endl;
01569   for (int k = 0 ; k < this->nbDomains; ++k) { // loop on subdomain
01570     //
01571     fileName = path + "/" + manager.getFiles().getResultFile() + "-density-"+ intToString(k) + ".ascii";
01572     //
01573     out.open(fileName.c_str());
01574     //
01575     out.setf(std::ios_base::scientific) ;   out.precision(20) ;
01576     //
01577     // ici il faut faire une boucle sur les SD
01578     //
01579     QCSymMatrix& QCRestrict  D = this->domains[k].getDensityP().getMatrix();
01580     //
01581     const int *globalIndex  = this->domains[k].getGlobalIdxs() ;
01582     QCSubDomainZone  * zone = this->domains[k].getZone();
01583     QCSubDomainZone  zonei, zonej ;
01584     QCPartitionType  partType = manager.getPartitioner()->getPartitionType() ;
01585     int cptA, cptB,flag ;
01586     //
01587     // Compute the number of elements to save
01588     //
01589     out  << "  "<< QCSystem::getNbAtomicOrbitals()  << "  "<< QCSystem::getNbAtomicOrbitals()  
01590          << "  "<<D.getNbElems() << std::endl; 
01591     //
01592     for (int i = 0; i <  this->domains[k].getNbAtoms(); ++i) {
01593     
01594       parameterA    = &parameters[ this->domains[k].getParamIndexAt(i)];
01595       firstAOofA    =  this->domains[k].getFirstAOAt(i);
01596       firstAOofGA   =  QCSystem::getFirstAOAt(globalIndex[i]); // accede au système global
01597       firstAOofSucA = firstAOofA + parameterA->getNbAO();
01598       zonei         = zone[i] ;
01599       //
01600       for (int j = 0; j <= i; ++j) {
01601         parameterB    = &parameters[ this->domains[k].getParamIndexAt(j)]; 
01602         firstAOofB    =  this->domains[k].getFirstAOAt(j); 
01603         firstAOofGB   =  QCSystem::getFirstAOAt(globalIndex[j]);
01604         firstAOofSucB = firstAOofB + parameterB->getNbAO();
01605         //
01606         zonej = zone[j] ;
01607         flag = computeFlag(zonei,zonej,partType); 
01608         cptA = firstAOofGA - 1 ;
01609         for (int mu = firstAOofA; mu < firstAOofSucA; ++mu) {
01610           ++cptA ;  cptB = firstAOofGB - 1 ;
01611           for (int nu = firstAOofB; nu < firstAOofSucB; ++nu) {
01612             ++cptB ;
01613             if( mu < nu) {  break;}
01614             if(cptA < cptB) {  //The global density matrix is also symtric
01615               out << "  "<< cptB << "  " << cptA <<  "  " << D[mu][nu] * flag << std::endl;
01616             }
01617             else {
01618               out << "  "<< cptA << "  " << cptB <<  "  " << D[mu][nu] *  flag<< std::endl;
01619             }
01620           }
01621         } 
01622       }
01623     }
01624     out.close();
01625   }
01626  // 
01627   QC_TRACE_OUT("END   QCMDSystem<TPSolver>::writeDensityOnFileAscii   " ) ;
01628 }
01629 template <class TPSolver>
01630 template <class TPManager>
01631 void 
01632 QCMDSystem<TPSolver>::writeOrbitalPartition(TPManager& manager,  const std::string& path){
01633   QC_TRACE_OUT("BEGIN QCMDSystem<TPSolver>::writeOrbitalPartition   ");
01634   //
01635   std::ofstream out;
01636   std::string fileName;
01637   //
01638   //  Il faut construire la correspondance entre l'orbitale locale et l'orbitale globale
01639   // Les parametres
01640   const typename TPManager::TParam *QCRestrict parameters =   manager.getParameters();
01641  // Le parametre courant
01642   const typename TPManager::TParam *QCRestrict parameterA;
01643   //
01644     //
01645   fileName = path + "/" + manager.getFiles().getResultFile() + "-orbital"+ ".part";
01646   //
01647   out.open(fileName.c_str());
01648   //
01649   QCPartitioner<QCMDSystem<TPSolver> > & QCRestrict  partitioner = (*manager.getPartitioner()) ;
01650   //
01651   out  << "Partitioning of the molecule "<< manager.getFiles().getAtomsFile() 
01652        << std::endl << std::endl ;
01653   partitioner.writeHeader(out);
01654   out << std::endl ;
01655   for (int k = 0 ; k < this->nbDomains; ++k) { // loop on subdomain
01656     const QCSubDomainZone  * zone         = this->domains[k].getZone() ;
01657     const int              * globalIndex  = this->domains[k].getGlobalIdxs() ;
01658     //
01659     int firstAOofA, firstAOofAG, nbAO, numCore, numShell1, numShell2  ;
01660     //
01661     numCore =  numShell1 = numShell2 = 0 ;
01662     for (int i = 0; i <  this->domains[k].getNbAtoms(); ++i) {
01663       parameterA    = &parameters[ this->domains[k].getParamIndexAt(i)];
01664       nbAO          = parameterA->getNbAO();
01665       //
01666       if (zone[i] == QC_CORE){
01667         numCore += nbAO; }
01668       else if (zone[i] == QC_SHELL1){
01669         numShell1 += nbAO; }
01670       else {
01671         numShell2 += nbAO ; }   
01672     }
01673     vector<int> core(numCore), shell1(numShell1),shell2(numShell2);
01674     vector<int> coreG(numCore), shell1G(numShell1),shell2G(numShell2);
01675     int ic,is1,is2 ;
01676     // Write orbitals in QC_SHELL1
01677     ic = is1 = is2 = 0;
01678     for (int i = 0; i <  this->domains[k].getNbAtoms(); ++i) {
01679       //  
01680       parameterA    = &parameters[ this->domains[k].getParamIndexAt(i)];
01681       //      firstAOofA    = this->domains[k].getFirstAOAt(i);
01682       nbAO        = parameterA->getNbAO();
01683       firstAOofAG  =  QCSystem::getFirstAOAt(globalIndex[i]);
01684       firstAOofA =  this->domains[k].getFirstAOAt(i);
01685       if (zone[i] == QC_CORE){
01686         for (int mu = 0; mu <  nbAO; ++mu) {
01687           core[ic]  = firstAOofA  + mu   ; 
01688           coreG[ic] = firstAOofAG + mu   ;++ic ;  
01689         }
01690       }   
01691       else if (zone[i] == QC_SHELL1){
01692         for (int mu = 0; mu <  nbAO; ++mu) {
01693           shell1[is1]  = firstAOofA  + mu   ;
01694           shell1G[is1] = firstAOofAG + mu   ; ++is1 ;  
01695         }
01696       }   
01697       else {
01698         for (int mu = 0; mu <  nbAO; ++mu) {
01699           shell2[is2]  = firstAOofA  + mu   ;
01700           shell2G[is2] = firstAOofAG + mu   ; ++is2 ;  
01701 
01702         }
01703       } 
01704     }
01705     // Write orbitals in QC_CORE
01706     out << "Sub Domain number : "<<k << " Number of Orbital :  "
01707         << this->domains[k].getNbAtomicOrbitals()  << std::endl ;
01708     out << "   CORE ZONE      : " <<core.size()<<std::endl <<"      ";
01709     for (unsigned int i = 0; i < core.size(); ++i) {
01710       //      out << core[i] << "  " ;}
01711       out << "( "<<core[i] <<", "<< coreG[i]<< ")  " ;}
01712     out <<std::endl << "   SHELL1 ZONE    : " <<shell1.size()<<std::endl << "   " ;
01713     for (unsigned int i = 0; i < shell1.size(); ++i) {
01714       out << "( "<<shell1[i] <<", "<< shell1G[i]<< ")  " ;}
01715     //      out << shell1[i] << "  " ;}
01716     out <<std::endl<< "   SHELL2 ZONE    : " <<shell2.size()<<std::endl<<  "   " ;
01717     for (unsigned int i = 0; i < shell2.size(); ++i) {
01718       out << "( "<<shell2[i]<< ", "<< shell2G[i]<< ")  " ;}
01719     //      out << shell2[i] << "  " ;}
01720     out <<std::endl<<std::endl;
01721   }
01722   out.close();
01723   // 
01724   QC_TRACE_OUT("END   QCMDSystem<TPSolver>::writeOrbitalPartition   " ) ;
01725 }
01726 template <class TPSolver>
01727 void 
01728 QCMDSystem<TPSolver>::writeSystemOnFile(const QCFiles & files,  const string& path){
01729   QC_TRACE_OUT("BEGIN QCMDSystem<TPSolver>::writeSystemOnFile   ") ;
01730   std::ofstream out ;
01731   std::string fileName;
01732   //
01733   fileName = path + "/" + files.getResultFile() + ".system";
01734   out.open(fileName.c_str());
01735   //
01736   QCSystem::writeQCIAtomicSystem(out) ;
01737 
01738   QC_TRACE_OUT("END   QCMDSystem<TPSolver>::writeLocalDensityOnFilesAscii  ") ;
01739 }
01740 
01741 
01742 //
01743 //
01744 template class QCMDSystem<QCDCAlgo>;
01745 
01746 template ostream& 
01747 operator << (ostream& out, const QCMDSystem<QCDCAlgo>& system);
01748 
01749 QCMANAGER_METH_EXPL_INST_MD_PARAM(void QCMDSystem<QCDCAlgo>::writeDensityOnFileAscii, ONE_PARAM(const std::string&));
01750 QCMANAGER_METH_EXPL_INST_MD_PARAM(void QCMDSystem<QCDCAlgo>::writeLocalDensityOnFilesAscii, ONE_PARAM(const std::string&));
01751 QCMANAGER_METH_EXPL_INST_MD_PARAM(void QCMDSystem<QCDCAlgo>::writeOrbitalPartition, ONE_PARAM(const std::string&));
01752 
01753 QCMANAGER_METH_EXPL_INST_MD_PARAM(void QCMDSystem<QCDCAlgo>::init, ONE_PARAM(const string&));
01754 QCMANAGER_METH_EXPL_INST_MD_PARAM(QCFloat QCMDSystem<QCDCAlgo>::completeFockMatricesAndElecEnergy, TWO_PARAMS(bool, int));
01755 
01756 
01757 QCMANAGER_METH_EXPL_INST_MD(void QCMDSystem<QCDCAlgo>::completeHamiltonMatrices);
01758 QCMANAGER_METH_EXPL_INST_MD(void QCMDSystem<QCDCAlgo>::completeDensityMatrices);
01759 
01760 QCPARAMETER_METH_EXPL_INST(QCFloat QCMDSystem<QCDCAlgo>::applyOptimalDamping);
01761 
01762 template void 
01763 QCMDSystem<QCDCAlgo>::buildPartitioner (QCPartitioner<TQCMDSystem>*&, const QCGeneralData&, const QCFiles&,
01764                                         const int *, const int);
01765 template void 
01766 QCMDSystem<QCDCAlgo>::allocatePartitioner (QCPartitioner<TQCMDSystem>*&, const QCGeneralData&); 
01767 
01768 template void
01769 QCMDSystem<QCDCAlgo>::adjustFermiEnergy(TQCMDSystem&, QCMemory&);
01770 
01771 template void
01772 QCMDSystem<QCDCAlgo>::adjustFermiEnergyNew(TQCMDSystem&, QCMemory&);
01773 
01774 
01775 #if defined(HAVE_MPI) && defined(WITH_MPI_SUPPORT)
01776 template void 
01777 QCMDSystem<QCDCAlgo>::allocatePartitioner (QCPartitioner<TQCDistMDSystem>*&, const QCGeneralData&); 
01778 template void 
01779 QCMDSystem<QCDCAlgo>::buildPartitioner (QCPartitioner<TQCDistMDSystem>*&, const QCGeneralData&,
01780                                         const QCFiles&, const int *,const int);
01781 template void
01782 QCMDSystem<QCDCAlgo>::adjustFermiEnergy(TQCDistMDSystem&, QCMemory&);
01783 
01784 template void
01785 QCMDSystem<QCDCAlgo>::adjustFermiEnergyNew(TQCDistMDSystem&, QCMemory&);
01786 //
01787 // 
01788 //
01789 QCMANAGER_METH_EXPL_INST_DIST_MD_PARAM(void QCMDSystem<QCDCAlgo>::writeDensityOnFileAscii, ONE_PARAM(const std::string&));
01790 
01791 #endif
01792 

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