QCDCAlgo.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 
00022 #include "QCMacro.hpp"
00023 #include "QCCommon.hpp"
00024 #include "QCManager.hpp"
00025 #include "QCReader.hpp"
00026 #include "QCMDSystem.hpp"
00027 
00028 #if defined(HAVE_MPI) && defined(WITH_MPI_SUPPORT) 
00029 #include "QCDistMDSystem.hpp"
00030 #   if defined(__QC_xlC__)
00031 #      include <VT_trc.h>
00032 #endif
00033 #endif
00034 #include "QCMndo.hpp"
00035 #include "QCMndoParam.hpp"
00036 #include "QCAm1.hpp"
00037 #include "QCAm1Param.hpp"
00038 #include "QCPm3.hpp"
00039 #include "QCPm3Param.hpp"
00040 #include "QCSCFAlgo.hpp"
00041 #include "QCChrono.hpp"
00042 //
00043 //
00044 #include "QCDCAlgo.hpp"
00045 //
00046 #if defined (QC_TRACE_IN_FILE)  
00047 #include "QCTrace.hpp"
00048 #endif
00049 
00050 
00054 QCDCAlgo::QCDCAlgo (void) : QCSCFBaseAlgo()
00055 {}
00059 QCDCAlgo::~QCDCAlgo (void)
00060 {}
00064 template <class TPManager>
00065 void QCDCAlgo::solve(TPManager& manager) {
00066   //
00067   std::cout.setf(ios::scientific);
00068   if (this->mainAlgo == QC_DC_SCF || mainAlgo == QC_DC_LEVEL_SHIFTING) {
00069     this->fixPoint (manager) ;
00070   }
00071   else if (this->mainAlgo == QC_DC_OPTIMAL_DAMPING) {
00072     this->optimalDamping (manager) ;
00073   } 
00074   else {
00075     std::cerr << __FILE__ << ": solve: method not implemented" << std::endl;
00076     exit(EXIT_FAILURE) ;
00077   }
00078 }
00082 template <class TPManager>
00083 void QCDCAlgo::fixPoint(TPManager& manager) {
00084   //
00085   QC_TRACE_SCF("BEGIN void QCDCAlgo::fixPoint(TPManager& manager)");
00086   //
00087   typedef typename TPManager::TSystem TSystem;
00088   typedef typename TPManager::TParam  TParam;
00089   //
00090   // Le gestionnaire de memoire
00091   //
00092   QCMemory& QCRestrict memory = manager.getMemory();
00093   //
00094   // Le systeme entier
00095   //
00096   TSystem& QCRestrict system = manager.getSystem();
00097   //
00098   //
00099   // Les parametres
00100   //
00101   const TParam * QCRestrict parameters = manager.getParameters();
00102   //
00103   // Un iterateur sur les domaines
00104   //
00105   typename TSystem::QCIterator current;
00106   //
00107   // l'identifiant du processus for the parallel version
00108   //
00109   int myrank  = system.getRank() ;
00110   bool master = (system.getRank() ==0 );
00111   //  
00112 #ifndef DEV_SCF_CLASS_NO_MEMORY
00113   //
00114   // densites de l'iteration precedente  P^{n-1} = Pnm1
00115   //
00116   QC_TRACE("            DEV_SCF_CLASS_NO_MEMORY est actif  ");
00117   for (current  = system.begin(); current != system.end();  current++) {
00118     current->getDensityPnm1().getMatrix().setDim(current->getDensityP().getMatrix().getDim());
00119     current->getDensityPnm1().getMatrix().takeTmpMem(memory);
00120   }
00121 #endif // DEV_SCF_CLASS_NO_MEMORY
00122   //
00123 #ifdef QC_VERBOSE_SCF
00124    if (master) { 
00125      std::cout << " ******************************************"  << std::endl
00126                << " ********** SOLVE: FIX POINT **************"  << std::endl
00127                << " ******************************************" << std::endl
00128                << std::endl;
00129    }
00130 #endif
00131   bool doLevelShifting = (this->mainAlgo == QC_DC_LEVEL_SHIFTING);
00132   QCFloat gap          = QC_ZERO;  
00133   QCFloat error        = 1000000 ;
00134   QCFloat elecEnergy, elecEnergyPrev;
00135 #if defined (QC_VERBOSE_SCF)
00136   if (master) {
00137     std::cout << std::endl << std::endl;
00138     std::cout << "PROC " << myrank << ":    **** [INIT/" << this->nbMaxIter << "] ****" << std::endl;
00139   }
00140 #endif
00141   //
00142   // Construction des matrices d'overlap et de Hamilton localement a chaque domaine
00143   //
00144   system.buildMatrices(manager);
00145   //
00146 #ifdef QC_REPORT_TIME
00147   // Les timer de mesures
00148   int i = 0;
00149   QCChrono   timer;
00150   QCFloat  * resolutionTimes = new QCFloat [system.getNbDomains()*nbMaxIter];
00151   QCFloat    tResTime        = 0.0;
00152   QCFloat  * densLTimes      = new QCFloat [nbMaxIter];
00153   QCFloat    tDensLTime      = 0.0;
00154   QCFloat  * fermiTimes      = new QCFloat [nbMaxIter];
00155   QCFloat    tFermiTime      = 0.0;
00156   QCFloat  * fockCTimes      = new QCFloat [nbMaxIter];
00157   QCFloat    tFockCTime      = 0.0;
00158   QCFloat  * fockLTimes      = new QCFloat [nbMaxIter];
00159   QCFloat    tFockLTime      = 0.0;
00160   QCFloat  * densCTimes      = new QCFloat [nbMaxIter];
00161   QCFloat    tDensCTime      = 0.0;
00162   ostringstream oss;
00163   oss << QCCommon::outdir << PATH_SEPARATOR << "TIME_DCSolver@" << myrank << ":" << QCCommon::getPidStr()
00164       << OUTPUT_SUFFIX;
00165   ofstream timelog (oss.str().c_str(), ios::out);  
00166 #endif
00167   //
00168   // Comms entre domaines pour la completion du hamiltonien pour chaque domaine
00169   //
00170   system.completeHamiltonMatrices(manager);
00171   //
00172 #ifdef QC_REPORT_TIME
00173   timelog << "*********** ITER 0 ***********"  << setprecision(10) << std::endl;
00174   timer.start();
00175 #endif
00176   //
00177   // On prepare les matrices de Fock en copiant les hamiltoniens
00178   //     et en completant localement a chaque domaine
00179   //
00180   system.prepareFockMatrices(manager, true);
00181   //
00182 #ifdef QC_REPORT_TIME
00183   timer.pause();
00184   fockLTimes[0] = timer.getvalsec();
00185   tFockLTime   += fockLTimes[0];
00186   timelog << " * fock local time:\t\t" << fockLTimes[0] << " s" << std::endl;
00187   timer.start();
00188 #endif
00189   //
00190   // Coms entre domaines pour la completion des matrices de fock
00191   //    et le calcul de la premiere energie
00192   //
00193   elecEnergy = system.completeFockMatricesAndElecEnergy(manager, true, 0);
00194   //
00195 #ifdef QC_REPORT_TIME
00196   timer.pause();
00197   fockCTimes[0] = timer.getvalsec();
00198   tFockCTime   += fockCTimes[0];
00199   timelog << " * fock contrib time:\t\t" << fockCTimes[0] << " s" << std::endl;
00200 #endif
00201   //  
00202   // On commence les Iterations
00203   //
00204   system.prepareIterations(manager);
00205   nbIter = 1;
00206   convReached = false;
00207   //
00208   // **************************************************************************
00209   // ****************** les iterations ****************************************
00210   // **************************************************************************
00211   //
00212   std::cout.setf(ios::scientific);
00213   while ( (!convReached) && (nbIter < nbMaxIter) ) {
00214     //
00215 #ifdef QC_REPORT_TIME      
00216     timelog << std::endl << "*********** ITER " << nbIter << " ***********" << std::endl;
00217 #endif
00218     //
00219     // Diagonalisation des matrices de fock et calcul du factorB
00220     //
00221     for (current  = system.begin(); current != system.end(); current++) {
00222       //
00223 #ifdef QC_REPORT_TIME      
00224       timer.start();
00225 #endif
00226       // Diagonalisation des matrices de fock
00227       current->getFockF().diagonalize(manager, current, nbIter);
00228       //
00229 #ifdef QC_REPORT_TIME
00230       timer.pause();
00231       resolutionTimes[i*nbMaxIter+nbIter]  = timer.getvalsec();
00232       tResTime                            += resolutionTimes[i*nbMaxIter+nbIter];
00233       timelog << " * fock diag time SD " << current->getId() << " :\t" << resolutionTimes[i*nbMaxIter+nbIter] 
00234               << " s" << std::endl;
00235       ++i;
00236 #endif
00237       //      
00238       // calcul du factorB
00239       current->getDensityPnm1().swapMatrix(current->getDensityP());
00240       current->computeBFactor(parameters);
00241     }
00242 #ifdef QC_REPORT_TIME
00243     timer.start();
00244 #endif
00245     //
00246     // On ajuste l'energie de fermi
00247     //
00248     //system.adjustFermiEnergy(system, memory);
00249     system.adjustFermiEnergyNew(system, memory);
00250     //
00251 #ifdef QC_REPORT_TIME
00252     timer.pause();
00253     fermiTimes[nbIter] = timer.getvalsec();
00254     tFermiTime        += fermiTimes[nbIter];
00255     timelog << " * fermi adjust time:\t\t" << fermiTimes[nbIter] << " s" << std::endl;
00256     i = 0;
00257     timer.start();
00258 #endif
00259     //
00260 #ifdef QC_DEBUG_SUMTRACE
00261     QCFloat sumTrace = 0;
00262 #endif
00263     //
00264     for (current  = system.begin(); current != system.end(); current++) {
00265       //
00266       current->getDensityP().computeElemsDC(memory, current->getNbFilledOM(), current->getNbOccupiedOM(),
00267                                             current->getNbAOCoreShell1(), current->getOrbitalOccupN(),
00268                                             current->getEigenVectCt());
00269       
00270       current->initWeightedDensity(parameters);
00271     }
00272 #ifdef QC_REPORT_TIME
00273     timer.pause();
00274     densLTimes[nbIter] = timer.getvalsec();
00275     tDensLTime        += densLTimes[nbIter];
00276     timelog << " * density local time:\t\t" << densLTimes[nbIter] << " s" << std::endl;
00277 #endif
00278 
00279     //
00280 #ifdef QC_DEBUG_SUMTRACE
00281     sumTrace = system.getTotalSumTrace();
00282     if (master) {
00283       std::cout << std::endl << "PROC " << myrank  << ": SUMTRACE = " << setprecision(10) << sumTrace << std::endl;
00284     }
00285 #endif
00286     //
00287 #ifdef QC_REPORT_TIME
00288     timer.start();
00289 #endif
00290     //
00291     // mise a jour de la matrice de densite avec les contributions des voisins
00292     //
00293     system.completeDensityMatrices(manager);
00294     //
00295 #ifdef QC_REPORT_TIME
00296     timer.pause();
00297     densCTimes[nbIter] = timer.getvalsec();
00298     tDensCTime        += densCTimes[nbIter];
00299     timelog << " * density contrib time:\t" << densCTimes[nbIter] << " s" << std::endl;
00300     timer.start();
00301 #endif
00302     //
00303     //  Build Fock matrix 
00304     //     Step 1 :  F = H
00305     //     Step 2 :  F = H + 1/2 G(D)
00306     //
00307     system.prepareFockMatrices(manager, false);
00308     //
00309 #ifdef QC_REPORT_TIME
00310     timer.pause();
00311     fockLTimes[nbIter] = timer.getvalsec();
00312     tFockLTime        += fockLTimes[nbIter];
00313     timelog << " * fock local time:\t\t" << fockLTimes[nbIter] << " s" << std::endl;
00314     timer.start();
00315 #endif
00316     //
00317     //  Contributions des voisins
00318     //
00319     //   La nouvelle energie 
00320     //
00321     elecEnergyPrev = elecEnergy;
00322     elecEnergy     = system.completeFockMatricesAndElecEnergy(manager, false, nbIter);
00323     //
00324 #ifdef QC_REPORT_TIME
00325     timer.pause();
00326     fockCTimes[nbIter] = timer.getvalsec();
00327     tFockCTime        += fockCTimes[nbIter];
00328     timelog << " * fock contrib time:\t\t" << fockCTimes[nbIter] << " s" << std::endl;
00329 #endif
00330     //
00331     //  Test de convergence 
00332     //
00333     convReached    = (system.testConvergence(error,threshold)&& (QCAbs(elecEnergy - elecEnergyPrev) <= threshold) );    
00334 #ifdef QC_VERBOSE_SCF
00335     if(master){
00336       std::cout << " ** iteration " << nbIter << "/" << nbMaxIter  
00337                 << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00338                 << std::setprecision(10) << elecEnergy  << std::endl ;
00339     }
00340 #endif
00341     //
00342     //
00343     // Le level shifting.
00344     //
00345     if (doLevelShifting) {
00346       std::cerr <<"  WARNING LEVEL SHIFTING NOT IMPLENETED IN DIVID AND CONQUER METHOD" <<std::endl;
00347       exit(1) ;
00348       //      system.computeLumoHomo(lumo,homo) ;
00349       //      gap  = homo - lumo;
00350       gap -= this->levelShiftingParam;
00351       //      fockF.levelShifting(workingSystem, gap);
00352 
00353     }
00354     //
00355     ++nbIter;
00356   }
00357   if(master){
00358     if(convReached) {
00359       std::cout <<std::endl << "  Convergence in "  ;}
00360     else {
00361       std::cout <<std::endl << "  No convergence in "  ;}
00362     std::cout << nbIter -1 << " iterations (" << nbMaxIter<<")"   
00363               << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00364               << std::setprecision(10) << elecEnergy  << std::endl << std::endl ;
00365   }
00366   
00367 #ifdef __QC_xlC__
00368   //VT_trc_stop_c ();    QCDensityGtr<QCSym    QCDensityGtr<QCSymMatrix>& QCRestrict Dtild = *matrices.getDensityPTild();
00369     QCFockGtr<QCSymMatrix>&    QCRestrict Ftild = *matrices.getFockFTild();
00370 Matrix>& QCRestrict Dtild = *matrices.getDensityPTild();
00371     QCFockGtr<QCSymMatrix>&    QCRestrict Ftild = *matrices.getFockFTild();
00372 
00373 #endif
00374 
00375 
00376 #ifdef QC_REPORT_TIME
00377   timelog << std::endl << std::endl
00378           << "#CUMULATED TIMES: Pid " << QCCommon::getPidStr()  << std::endl << std::endl 
00379           << " * fock diag time:\t\t"     << tResTime << " s"   << std::endl 
00380           << " * fermi adjust time:\t\t"  << tFermiTime << " s" << std::endl  
00381           << " * density local time:\t\t" << tDensLTime << " s" << std::endl
00382           << " * density contrib time:\t" << tDensCTime << " s" << std::endl
00383           << " * fock local time:\t\t"    << tFockLTime << " s" << std::endl
00384           << " * fock contrib time:\t\t"  << tFockCTime << " s" << std::endl; 
00385 
00386    system.computeAverages(tResTime,
00387                          tFermiTime,
00388                          tDensLTime,
00389                          tDensCTime,
00390                          tFockLTime,
00391                          tFockCTime);
00392 
00393    if (system.getRank() == 0) {
00394      timelog << std::endl << std::endl   << "# AVERAGE TIMES:"  << std::endl << std::endl
00395              << " * fock diag time (per domain):\t\t"  << tResTime << " s" << std::endl
00396              << " * fermi adjust time:\t\t\t"   << tFermiTime << " s" << std::endl
00397             << " * density local time (per domain):\t"    << tDensLTime << " s" << std::endl
00398             << " * density contrib time:\t\t"    << tDensCTime << " s" << std::endl
00399             << " * fock local time (per domain):\t"     << tFockLTime << " s" << std::endl
00400             << " * fock contrib time:\t\t\t"    << tFockCTime << " s" << std::endl;
00401    }
00402 
00403   timelog.close();
00404   //
00405   delete [] resolutionTimes;
00406   delete [] densLTimes;
00407   delete [] densCTimes;
00408   delete [] fermiTimes;
00409   delete [] fockCTimes;
00410   delete [] fockLTimes;
00411 #endif
00412 
00413 
00414 #ifndef DEV_SCF_CLASS_NO_MEMORY
00415   // On libere l'espace memoire pour les matrices
00416   // densites de l'iteration precedente
00417   for (current  = system.end();current != system.begin(); current--) {
00418      current->getDensityPnm1().getMatrix().giveBackTmpMem(memory);
00419   }
00420 #endif // DEV_SCF_CLASS_NO_MEMORY
00421   QC_TRACE("END void QCDCAlgo::fixPoint(TPManager& manager)");
00422 }
00426 #if defined( BUILD_OC)
00427 template <class TPManager>
00428 void QCDCAlgo::calcHomoLumo (const TPManager& manager, QCFloat & globalHomo, QCFloat & globalLumo) {
00429   //
00430   typedef typename TPManager::TSystem TSystem;
00431   //
00432   TSystem& QCRestrict system = manager.getSystem();
00433   typename TSystem::QCIterator current;
00434   //
00435   const DBL_MAX = 1.0e23 ;
00436   //
00437   // L'homo et la lumo du SD courant
00438   //
00439   QCFloat homoOfSD, lumoOfSD;
00440   //
00441   int nbOccupiedOM, numSD;
00442   //
00443   //La plus haute des homo des sous systemes.
00444   //
00445   globalHomo = -DBL_MAX;
00446   globalLumo = DBL_MAX;
00447   //
00448   for (numSD = 0; numSD < nbSubDomains; ++numSD) {
00449     nbOccupiedOM = subDomains[numSD].getSubModel().getNbOccupiedOM();
00450     //
00451     //  Attention, on verifie qu il existe bien une homo et que le sous-domaine
00452     //      n'est pas completement vide.
00453     //
00454     homoOfSD = -DBL_MAX;
00455     if (nbOccupiedOM > 0) {
00456       homoOfSD = eigenValArray[numSD][nbOccupiedOM-1];
00457     }
00458     if (homoOfSD > globalHomo) {
00459       globalHomo = homoOfSD;
00460     }
00461     //
00462     // Attention, on verifie qu'il existe bien une Lumo et que le sous-domaine
00463     //    n'est pas completement rempli.
00464     //
00465     lumoOfSD = DBL_MAX;
00466     if (nbOccupiedOM < subDomains[numSD].getSubSystem().getNbOMOfSystem() ) {
00467       lumoOfSD = eigenValArray[numSD][nbOccupiedOM];
00468     }
00469     if (lumoOfSD < globalLumo) {
00470       globalLumo = lumoOfSD;
00471     }
00472   }//for numSD
00473 
00474 }
00475 template <class TPManager>
00476 void QCDCAlgo::applyLevelShifting (TPManager& managerInst) {
00477   //
00478   //  Le gap entre la lumo et l homo du systeme global.
00479   QCFloat globalGap, globalLumo, globalHomo ;
00480   //
00481   int numSD;
00482   //
00483   //Calcul de l homo et de la lumo.
00484   this->calcHomoLumo(globalLumo, globalHomo);
00485   //
00486   globalGap = globalLumo - globalHomo;
00487   cout << "globalGap = " << globalGap << endl;
00488   //
00489   //Le shift de level si necessaire.
00490   //
00491   if (globalGap < this->levelShiftingParam ) {
00492     for (numSD = 0; numSD < nbSubDomains; ++numSD) {
00493       managerInst.setComputedData(subDomains[numSD]);
00494       subDomains[numSD].getSubModel().getFockF().levelShifting(managerInst, globalGap);
00495     }//for numSD
00496   }//if gap < managerInst.getShiftingLevelParameter
00497 }
00498 
00499 #endif
00500 
00501 template <class TPManager>
00502 void QCDCAlgo::optimalDamping(TPManager& manager) {
00503   //
00504   QC_TRACE_SCF("BEGIN void QCDCAlgo::optimalDamping(TPManager& manager)");
00505   //
00506   typedef typename TPManager::TSystem TSystem;
00507   typedef typename TPManager::TParam  TParam;
00508   //
00509   //
00510   // Le gestionnaire de memoire
00511   //
00512   QCMemory& QCRestrict memory = manager.getMemory();
00513   //
00514   // Le systeme entier
00515   //
00516   TSystem& QCRestrict system = manager.getSystem();
00517   //
00518   // Les parametres
00519   //
00520   const TParam * QCRestrict parameters = manager.getParameters();
00521   //
00522   // Un iterateur sur les domaines
00523   //
00524   typename TSystem::QCIterator current;
00525   //
00526   // l'identifiant du processus for the parallel version
00527   //
00528   int  myrank = system.getRank() ;
00529   bool master = (myrank ==0 );
00530   //
00531   // Les matrices temporaires
00532   //
00533 #ifndef DEV_SCF_CLASS_NO_MEMORY
00534   // On aloue l'espace memoire pour les matrices via le gestionnaire de memoire 
00535   //      - densites de l'iteration precedente Pn-1
00536   //      - Ftild et Dtild
00537   int dim ;
00538   for (current  = system.begin(); current != system.end();  current++) { 
00539     dim = current->getDensityP().getMatrix().getDim() ;
00540     current->getDensityPnm1().getMatrix().setDim(dim); current->getDensityPnm1().getMatrix().takeTmpMem(memory);
00541     current->getDensityPTild()->getMatrix().setDim(dim); current->getDensityPTild()->.getMatrix().takeTmpMem(memory);
00542     current->getFockFTild()->getMatrix().setDim(dim); current->getFockFTild()->getMatrix().takeTmpMem(memory);
00543   }
00544 #endif // DEV_SCF_CLASS_NO_MEMORY
00545   //
00546 #ifdef  QC_VERBOSE_GNUPLOT 
00547   std::ofstream gnuplot("scf.plot") ;
00548 #endif
00549 #ifdef QC_VERBOSE_SCF
00550   if (master) {
00551     std::cout << " ************************************************"  << std::endl
00552               << " ********** SOLVE: Optimal Damping **************"  << std::endl
00553               << " ************************************************" << std::endl
00554               << std::endl;
00555     std::cout << std::endl << std::endl;
00556     std::cout << "PROC " << myrank << ":    **** [INIT/" << this->nbMaxIter << "] ****" << std::endl;
00557   }
00558 #endif
00559 #if defined (QC_TRACE_IN_FILE)  
00560   int ii ;
00561   for (current  = system.begin(), ii=0; current != system.end();  ++current,++ii  ) { 
00562     qctrace[ii].writeMsg("Optimal Damping Solver");
00563   }
00564 #endif
00565 #ifdef QC_REPORT_TIME
00566   // Les timer de mesures
00567   int i = 0;
00568   QCChrono   timer;
00569   QCFloat  * resolutionTimes = new QCFloat [system.getNbDomains()*nbMaxIter];
00570   QCFloat    tResTime        = 0.0;
00571   QCFloat  * densLTimes      = new QCFloat [nbMaxIter];
00572   QCFloat    tDensLTime      = 0.0;
00573   QCFloat  * fermiTimes      = new QCFloat [nbMaxIter];
00574   QCFloat    tFermiTime      = 0.0;
00575   QCFloat  * fockCTimes      = new QCFloat [nbMaxIter];
00576   QCFloat    tFockCTime      = 0.0;
00577   QCFloat  * fockLTimes      = new QCFloat [nbMaxIter];
00578   QCFloat    tFockLTime      = 0.0;
00579   QCFloat  * densCTimes      = new QCFloat [nbMaxIter];
00580   QCFloat    tDensCTime      = 0.0;
00581   ostringstream oss;
00582   oss << QCCommon::outdir << PATH_SEPARATOR << "TIME_DCSolver@" << myrank << ":" << QCCommon::getPidStr()
00583       << OUTPUT_SUFFIX;
00584   ofstream timelog (oss.str().c_str(), ios::out);  
00585 #endif
00586   //
00587   // Les energies
00588   //
00589   //  QCFloat coreEnergy, elecEnergy, elecEnergyPrev, totalEnergy;
00590   QCFloat elecEnergy , elecEnergyPrev ;//, energy_tild, energy_tild_1e, energy_1e ; //, totalEnergy;
00591 #ifdef DEV_OL
00592   QCFloat c1,t, s, c, lambdaOpt ;
00593   QCFloat energy_tild_2e, energy_tild, energy_tild_1e, energy_1e;
00594 #endif
00595   elecEnergy = elecEnergyPrev = QC_ZERO ; 
00596 //   //
00597 //   int nbElectrons   = workingSystem->getNbElectrons();    
00598 //  int nbOccupiedOM  ;
00599 //
00600    // Construction des matrices d'overlap et de Hamilton localement a chaque domaine
00601   //
00602   system.buildMatrices(manager);
00603   //
00604   // Comms entre domaines pour la completion du hamiltonien pour chaque domaine
00605   //
00606   system.completeHamiltonMatrices(manager);
00607   
00608 #if defined (QC_TRACE_IN_FILE)  
00609   for (current  = system.begin(),ii=0; current != system.end();  ++current,++ii  ) { 
00610     QCModelMatrices& QCRestrict matrices = current->getMatrices();
00611     qctrace[ii].writeHamiltonMatrix(matrices);
00612     qctrace[ii].writeDensityMatrix(matrices);
00613     qctrace[ii].writeMsg(" ++++++++ END initialization ++++");
00614   }
00615 #endif
00616   //
00617 #ifdef QC_REPORT_TIME
00618   timelog << "*********** ITER 0 ***********"  << setprecision(10) << std::endl;
00619   timer.start();
00620 #endif
00621   //
00622   // On prepare les matrices de Fock en copiant les hamiltoniens
00623   //     et en completant localement a chaque domaine
00624   //
00625   system.prepareFockMatrices(manager, true);
00626   //
00627 #ifdef QC_REPORT_TIME
00628   timer.pause();
00629   fockLTimes[0] = timer.getvalsec();
00630   tFockLTime   += fockLTimes[0];
00631   timelog << " * fock local time:\t\t" << fockLTimes[0] << " s" << std::endl;
00632   timer.start();
00633 #endif
00634   //
00635   // Coms entre domaines pour la completion des matrices de fock
00636   //    et le calcul de la premiere energie
00637   //
00638   elecEnergyPrev = elecEnergy ;
00639   elecEnergy     = system.completeFockMatricesAndElecEnergy(manager, true, 0);
00640   //
00641 #ifdef QC_REPORT_TIME
00642   timer.pause();
00643   fockCTimes[0] = timer.getvalsec();
00644   tFockCTime   += fockCTimes[0];
00645   timelog << " * fock contrib time:\t\t" << fockCTimes[0] << " s" << std::endl;
00646 #endif
00647   //
00648   system.prepareIterations(manager);
00649   //
00650   // Calcul des energies du systeme  
00651   //
00652   //   coreEnergy  = model.computeCoreEnergy (manager, workingSystem);
00653   //elecEnergy  = model.computeElecEnergy (workingSystem, densityP);
00654   //   totalEnergy = coreEnergy + elecEnergy; 
00655   //   workingSystem->setEnergy (totalEnergy);
00656 
00657   nbIter      = 0;
00658   convReached = false;
00659   QCFloat error;
00660   //
00661   //   t = s = c = lambdaOpt = QC_ZERO;
00662   //   energy_1e      = QC_HALF * H.traceProduct(D) ;
00663   //  energy_tild    = elecEnergy;
00664   //   energy_tild_1e = energy_1e;
00665   //
00666 #if defined (QC_TRACE_IN_FILE)  
00667   for (current  = system.begin(), ii=0; current != system.end();  ++current ,++ii ) { 
00668 #else
00669   for (current  = system.begin(); current != system.end();  ++current  ) { 
00670 #endif
00671     QCModelMatrices& QCRestrict matrices = current->getMatrices();
00672     //
00673     // On initialise F_tilde avec F
00674     //
00675     matrices.getFockFTild()->copyMatrix(matrices.getFockF());
00676 #if defined (QC_TRACE_IN_FILE)  
00677     qctrace[ii].writeFockMatrix(matrices);
00678 #endif
00679     //
00680     // On initialise P_tilde avec P
00681     //
00682     matrices.getDensityPTild()->copyMatrix(matrices.getDensityP());
00683   }
00684   //
00685   nbIter = 1;
00686   convReached = false;
00687   //
00688   // **************************************************************************
00689   // ****************** les iterations ****************************************
00690   // **************************************************************************
00691   //
00692   while ( (!convReached) &&  (nbIter < nbMaxIter) ) {
00693     //    
00694 #ifdef QC_REPORT_TIME      
00695     timelog << std::endl << "*********** ITER " << nbIter << " ***********" << std::endl;
00696     timer.start();
00697 #endif
00698 #if defined (QC_TRACE_IN_FILE)  
00699   for (current  = system.begin(),  ii=0; current != system.end();  ++current,++ii  ) { 
00700     if(nbIter > 2 && nbIter < 5) {qctrace[ii].traceOn();}  else   {qctrace[ii].traceOff();}
00701     qctrace[ii].writeVal("  Iteration number : ",nbIter);
00702   }
00703   ii = 0;
00704 #endif
00705     //
00706     // Diagonalisation des matrices de focktild et calcul du factorB
00707     //
00708   for (current  = system.begin(); current != system.end(); current++) {
00709       //
00710 #ifdef QC_REPORT_TIME      
00711       timer.start();
00712 #endif
00713       // Diagonalisation des matrices de fock
00714       current->getFockFTild()->diagonalize(manager, current, nbIter);
00715       //
00716 #ifdef QC_REPORT_TIME
00717       timer.pause();
00718       resolutionTimes[i*nbMaxIter+nbIter]  = timer.getvalsec();
00719       tResTime                            += resolutionTimes[i*nbMaxIter+nbIter];
00720       timelog << " * fock diag time SD " << current->getId() << " :\t" << resolutionTimes[i*nbMaxIter+nbIter] 
00721               << " s" << std::endl;
00722       ++i;
00723 #endif
00724       //      
00725       // calcul du factorB
00726       current->getDensityPnm1().swapMatrix(current->getDensityP()); 
00727       current->computeBFactor(parameters);
00728 #if defined (QC_TRACE_IN_FILE)  
00729       //      qctrace[ii].writeMsg("Density Matrix P at time n-1") ;
00730       qctrace[ii].writeEigenValuesAndEigenVectors();
00731       ++ii;
00732 #endif
00733 
00734     }
00735   elecEnergyPrev = elecEnergy;
00736 #ifdef QC_REPORT_TIME
00737     timer.start();
00738 #endif
00739     //
00740     // On ajuste l'energie de fermi
00741     //
00742     //    system.adjustFermiEnergy(system, memory);
00743     system.adjustFermiEnergySeq(system, memory);
00744     //    system.adjustFermiEnergyNew(system, memory);
00745     //
00746 #ifdef QC_REPORT_TIME
00747     timer.pause();
00748     fermiTimes[nbIter] = timer.getvalsec();
00749     tFermiTime        += fermiTimes[nbIter];
00750     timelog << " * fermi adjust time:\t\t" << fermiTimes[nbIter] << " s" << std::endl;
00751     i = 0;
00752     timer.start();
00753 #endif
00754     //
00755 #ifdef QC_DEBUG_SUMTRACE
00756     QCFloat sumTrace = 0;
00757 #endif
00758     //
00759 #if !defined (QC_TRACE_IN_FILE)  
00760     for (current  = system.begin(); current != system.end(); ++current) {
00761 #else
00762     for (current  = system.begin(), ii=0; current != system.end(); ++current,++ii) {
00763       qctrace[ii].writeVal("Nb Filled OM        : ",current->getNbFilledOM());
00764       qctrace[ii].writeVal("Nb Occupied OM      : ",current->getNbOccupiedOM());
00765       qctrace[ii].writeVal("Nb AO Core Shell1   : ",current->getNbAOCoreShell1());
00766 //       qctrace[ii].writeArray(current->getNbFilledOM(),current->getNbOccupiedOM(),
00767 //                           current->getOrbitalOccupN());
00768 #endif
00769       current->getDensityP().computeElemsDC(memory, current->getNbFilledOM(), current->getNbOccupiedOM(),
00770                                             current->getNbAOCoreShell1(), current->getOrbitalOccupN(),
00771                                             current->getEigenVectCt());
00772       
00773       current->initWeightedDensity(parameters);
00774 #if defined (QC_TRACE_IN_FILE)  
00775       qctrace[ii].writeDensityMatrix();
00776       qctrace[ii].writeReducedDensityMatrix();
00777 #endif
00778     }
00779 #ifdef QC_REPORT_TIME
00780     timer.pause();
00781     densLTimes[nbIter] = timer.getvalsec();
00782     tDensLTime        += densLTimes[nbIter];
00783     timelog << " * density local time:\t\t" << densLTimes[nbIter] << " s" << std::endl;
00784 #endif
00785 
00786     //
00787 #ifdef QC_DEBUG_SUMTRACE
00788     sumTrace = system.getTotalSumTrace();
00789     if (master) {
00790       std::cout << std::endl << "PROC " << myrank  << ": SUMTRACE = " << setprecision(10) << sumTrace << std::endl;
00791     }
00792 #endif
00793     //
00794 #ifdef QC_REPORT_TIME
00795     timer.start();
00796 #endif
00797     //
00798     // mise a jour de la matrice de densite avec les contributions des voisins
00799     //
00800     system.completeDensityMatrices(manager);
00801     //
00802 #ifdef QC_REPORT_TIME
00803     timer.pause();
00804     densCTimes[nbIter] = timer.getvalsec();
00805     tDensCTime        += densCTimes[nbIter];
00806     timelog << " * density contrib time:\t" << densCTimes[nbIter] << " s" << std::endl;
00807     timer.start();
00808 #endif
00809     //    
00810     //  Build Fock matrix 
00811     //     Step 1 :  F = H
00812     //     Step 2 :  F = H + 1/2 G(D)
00813     //
00814     system.prepareFockMatrices(manager, false);
00815     //
00816 #if defined (QC_TRACE_IN_FILE)  
00817     for (current  = system.begin(), ii=0; current != system.end(); ++current,++ii) {
00818       qctrace[ii].writeMsg("mise a jour de la matrice de densite avec les contributions des voisins");
00819       qctrace[ii].writeDensityMatrix();
00820       qctrace[ii].writeMsg("Local Fock Matrix");
00821       qctrace[ii].writeFockMatrix() ;
00822       qctrace[ii].writeMsg("Add contribution of all subdomains -  Fock Matrix");
00823     }
00824 #endif
00825 #ifdef QC_REPORT_TIME
00826     timer.pause();
00827     fockLTimes[nbIter] = timer.getvalsec();
00828     tFockLTime        += fockLTimes[nbIter];
00829     timelog << " * fock local time:\t\t" << fockLTimes[nbIter] << " s" << std::endl;
00830     timer.start();
00831 #endif
00832     //
00833     //  Contributions des voisins
00834     //
00835     elecEnergy = system.completeFockMatricesAndElecEnergy(manager, false, nbIter);
00836     //
00837     //  energy_1e = QC_HALF * H.traceProduct(D);   A Faire !!!
00838 #ifdef QC_REPORT_TIME
00839     timer.pause();
00840     fockCTimes[nbIter] = timer.getvalsec();
00841     tFockCTime        += fockCTimes[nbIter];
00842     timelog << " * fock contrib time:\t\t" << fockCTimes[nbIter] << " s" << std::endl;
00843 #endif
00844     //
00845     // *******************************
00846     // ** Minimisation de l'energie **
00847     // *******************************
00848     //
00849 #ifdef DEV_OL
00850     lambdaOpt = QC_ONE ;
00851     t = system.computeSumOfEigenValOccupied() ;
00852     energy_tild_2e  = energy_tild - QC_TWO * energy_tild_1e ;
00853     s               = QC_TWO * t - QC_TWO * (energy_tild_1e + energy_tild_2e);
00854     c               = elecEnergy + energy_tild_2e - QC_TWO * t;
00855     //
00856     // lambdaOpt = min { -s / QCTwo * c), 1}
00857     //
00858     if (c > QC_ZERO && -s < (QC_TWO * c) ) {
00859       lambdaOpt = -s / (QC_TWO * c);
00860     }
00861     c1 = QC_ONE - lambdaOpt ;
00862     std::cout <<" t : " << t << " lambda Opt " << lambdaOpt << " -s / QCTwo * c "<< -s / QC_TWO * c<< std::endl;
00863     //
00864     //   update of Ftild and Dtild
00865     //
00866     //     elecEnergyPrev = elecEnergy;
00867     //     elecEnergy     = model.computeElecEnergy(workingSystem, densityP);
00868     //     totalEnergy    = coreEnergy + elecEnergy;
00869     //     energy_1e      = QC_HALF * H.traceProduct(D);
00870     //
00871     //
00872     // energy_tild update 
00873     //
00874     energy_tild    += lambdaOpt * ( s + lambdaOpt *c) ;
00875     energy_tild_1e  = c1 * energy_tild_1e + lambdaOpt * energy_1e ;
00876 #endif
00877     //
00878     //  Test de convergence 
00879     //
00880     convReached = (system.testConvergence(error,threshold) &&(QCAbs(elecEnergy-elecEnergyPrev)<=threshold));
00881     //
00882 #ifdef QC_VERBOSE_SCF
00883     if(master){
00884       std::cout << " ** iteration " << nbIter << "/" << nbMaxIter  
00885                 << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00886                 << std::setprecision(10) << elecEnergy  << std::endl ;
00887     }
00888 #endif
00889     //
00890     //  Minimisation de l'energie 
00891     //  Update : matrices :Ftild, Dtild and energy_tild 
00892     //
00893 #ifdef  QC_VERBOSE_GNUPLOT 
00894     QCFloat lambdaOpt =
00895 #endif
00896       system.applyOptimalDamping(parameters) ;
00897     //
00898 #ifdef  QC_VERBOSE_GNUPLOT 
00899     if(master){
00900       gnuplot << nbIter << "  "<< log(error) << " " << lambdaOpt <<std::endl;
00901     }
00902 #endif
00903     ++nbIter;
00904     }
00905     if(master){
00906       if(convReached) {
00907         std::cout <<std::endl << "  Convergence in "  ;}
00908       else {
00909         std::cout <<std::endl << "  No convergence in "  ;}
00910       std::cout << nbIter -1 << " iterations (" << nbMaxIter<<")"   
00911                 << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00912                 << std::setprecision(10) << elecEnergy  << std::endl << std::endl ;
00913     }
00914   //
00915   //  La matrice de densite est-elle indempotente ?
00916   //
00917   //   QCFloat indempo = D.traceProduct(D) - QCTWO * D.trace() ;
00918   //   if (convReached  && QCAbs(indempo) > threshold) {
00919   //     convReached = false;
00920   //     std::cout << "\n WARNING: The demsity matrix is not indempotent !!\n" << std::endl;
00921   //   }
00922 #ifdef QC_REPORT_TIME
00923   timelog << std::endl << std::endl
00924           << "#CUMULATED TIMES: Pid " << QCCommon::getPidStr()  << std::endl << std::endl 
00925           << " * fock diag time:\t\t"     << tResTime << " s"   << std::endl 
00926           << " * fermi adjust time:\t\t"  << tFermiTime << " s" << std::endl  
00927           << " * density local time:\t\t" << tDensLTime << " s" << std::endl
00928           << " * density contrib time:\t" << tDensCTime << " s" << std::endl
00929           << " * fock local time:\t\t"    << tFockLTime << " s" << std::endl
00930           << " * fock contrib time:\t\t"  << tFockCTime << " s" << std::endl; 
00931   timelog.close();
00932   //
00933   delete [] resolutionTimes;
00934   delete [] densLTimes;
00935   delete [] densCTimes;
00936   delete [] fermiTimes;
00937   delete [] fockCTimes;
00938   delete [] fockLTimes;
00939   //
00940 #endif
00941   //
00942 #ifndef DEV_SCF_CLASS_NO_MEMORY
00943   // On liberel'espace memoire du gestionnaire de memoire 
00944   //      - densites de l'iteration precedente Pn-1
00945   //      - Ftild et Dtild
00946   for (current  = system.begin(); current != system.end();  current++) { 
00947     current->getDensityPnm1().getMatrix().giveBackTmpMem(memory);
00948     current->getDensityPTild()->getMatrix().giveBackTmpMem(memory);
00949     current->getFockFTild()->getMatrix().giveBackTmpMem(memory);
00950   }
00951 #endif // DEV_SCF_CLASS_NO_MEMORY
00952 
00953   QC_TRACE("END   void QCDCAlgo::optimalDamping(TPManager& manager)");
00954 }
00955 
00956 
00957 QCMANAGER_METH_EXPL_INST_MD(void QCDCAlgo::solve);
00958 QCMANAGER_METH_EXPL_INST_MD(void QCDCAlgo::fixPoint);
00959 QCMANAGER_METH_EXPL_INST_MD(void QCDCAlgo::optimalDamping);
00960 
00961 #if defined(HAVE_MPI) && defined(WITH_MPI_SUPPORT)
00962 
00963 QCMANAGER_METH_EXPL_INST_DIST_MD(void QCDCAlgo::solve);
00964 QCMANAGER_METH_EXPL_INST_DIST_MD(void QCDCAlgo::fixPoint);
00965 QCMANAGER_METH_EXPL_INST_DIST_MD(void QCDCAlgo::optimalDamping);
00966 #endif

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