QCSCFAlgo.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 #include <iomanip>
00020 
00021 #include "QCMacro.hpp"
00022 #include "QCCommon.hpp"
00023 #include "QCManager.hpp"
00024 #include "QCReader.hpp"
00025 #include "QCGlobalSystem.hpp"
00026 #include "QCMDSystem.hpp"
00027 #include "QCMndo.hpp"
00028 #include "QCMndoParam.hpp"
00029 #include "QCAm1.hpp"
00030 #include "QCAm1Param.hpp"
00031 #include "QCPm3.hpp"
00032 #include "QCPm3Param.hpp"
00033 #include "QCSCFAlgo.hpp"
00034 #include "QCChrono.hpp"
00035 
00036 #if defined (QC_TRACE_IN_FILE)  
00037 #include "QCTrace.hpp"
00038 #endif
00039 
00042 QCSCFAlgo::QCSCFAlgo (void) : QCSCFBaseAlgo()
00043 {}
00044 
00048 QCSCFAlgo::~QCSCFAlgo (void)
00049 {}
00053 void QCSCFAlgo::setTempMemoryCost (const int /* dim */) {
00054 }
00055 //
00056 //
00057 //
00058 template <class TPManager>
00059 void QCSCFAlgo::solve (TPManager& manager) {
00060   //  
00061   std::cout.setf(ios::scientific);
00062   if (mainAlgo == QC_SCF || mainAlgo == QC_LEVEL_SHIFTING) {
00063     fixPoint (manager) ;
00064   }
00065   else if (mainAlgo == QC_OPTIMAL_DAMPING) {
00066     optimalDamping (manager) ;
00067   } 
00068   else {
00069     std::cerr << __FILE__ << ": solve: method not implemented" << std::endl;
00070     exit(EXIT_FAILURE) ;
00071   }
00072 
00073 }
00074 
00078 template <class TPManager>
00079 void QCSCFAlgo::fixPoint (TPManager& manager) {
00080 
00081 #ifdef QC_VERBOSE_SCF
00082   std::cout << " ******************************************"  << std::endl
00083             << " ********** SOLVE: FIX POINT **************"  << std::endl
00084             << " ******************************************" << std::endl
00085             << std::endl;
00086 #endif
00087 
00088 #ifdef QC_REPORT_TIME
00089   // Les timer de mesures
00090   QCChrono   timer;
00091   QCFloat  * diagTimes       = new QCFloat [nbMaxIter];
00092   QCFloat    tDiagTime       = 0.0;
00093   QCFloat  * densTimes       = new QCFloat [nbMaxIter];
00094   QCFloat    tDensTime       = 0.0;
00095   QCFloat  * fockTimes       = new QCFloat [nbMaxIter];
00096   QCFloat    tFockTime       = 0.0;
00097   string logname;
00098   logname = ( QCCommon::outdir + PATH_SEPARATOR + "TIME_SCFSolver@" + QCCommon::getPidStr()
00099               + OUTPUT_SUFFIX );
00100   ofstream timelog (logname.c_str(), ios::out);  
00101 #endif
00102   //
00103   // le gestionnaire de memoire
00104   //
00105   QCMemory& QCRestrict memory = manager.getMemory();
00106   //
00107   // Le systeme entier
00108   //
00109   typename TPManager::TSystem& QCRestrict system =  manager.getSystem();
00110   //
00111   // Le modele
00112   //
00113   typename TPManager::TModel& model = manager.getModel();
00114   //
00115   // L'iterateur
00116   //
00117   typename TPManager::TSystem::QCIterator workingSystem =  system.begin();
00118   //  
00119   // les matrices courantes
00120   //
00121   QCModelMatrices& QCRestrict matrices = workingSystem->getMatrices();
00122   //
00123   // Les alias vers les matrices
00124   //
00125   QCFockGtr<QCSymMatrix>&    QCRestrict fockF     = matrices.getFockF();
00126   QCDensityGtr<QCSymMatrix>& QCRestrict densityP  = matrices.getDensityP();
00127   QCMatElemGtr<QCMatrix>&    QCRestrict eigenVect = matrices.getEigenVectCt();
00128   //
00129   QCSymMatrix& QCRestrict  D        = densityP.getMatrix();
00130   QCSymMatrix& QCRestrict  H        = matrices.getHamiltonH().getMatrix();
00131   QCFloat *QCRestrict      eigenVal = matrices.getEigenVal();
00132   //
00133   // Les energies
00134   //
00135   QCFloat coreEnergy, elecEnergy, elecEnergyPrev, totalEnergy;
00136   //
00137   // La matrice de densite precedente
00138   //
00139   QCSymMatrix& QCRestrict Dprev = matrices.getDensityPnm1().getMatrix(); 
00140 
00141 #ifndef DEV_SCF_CLASS_NO_MEMORY
00142   QC_TRACE("            DEV_SCF_CLASS_NO_MEMORY est actif  ");
00143   Dprev.setDim(D.getDim());    Dprev.takeTmpMem(memory); 
00144 #endif
00145   //
00146   bool doLevelShifting = (this->mainAlgo == QC_LEVEL_SHIFTING);
00147   QCFloat error        = 1000000 ;
00148   QCFloat gap ;
00149 
00150   int nbElectrons      = workingSystem->getNbElectrons();    
00151   int nbOccupiedOM     = nbElectrons / 2; 
00152   //
00153   // Construction des matrices
00154   //
00155   system.buildMatrices(manager);
00156   //
00157   // Calcul des energies du systeme  
00158   //
00159   coreEnergy  = model.computeCoreEnergy (manager, workingSystem);
00160   elecEnergy  = model.computeElecEnergy (workingSystem, densityP);
00161   totalEnergy = coreEnergy + elecEnergy; 
00162 
00163 #ifdef QC_VERBOSE_SCF
00164   std::cout << setprecision(20)  << std::endl << "  * Electronic energy = " << elecEnergy  << std::endl
00165             << "  * Total energy      = " << totalEnergy << std::endl;
00166 #endif
00167   //
00168   workingSystem->setEnergy (totalEnergy);
00169   //
00170   // On commence les Iterations
00171   //
00172   system.prepareIterations(manager);
00173   nbIter = 0;
00174   convReached = false;
00175   //
00176   // **************************************************************************
00177   // ****************** les iterations ****************************************
00178   // **************************************************************************
00179   //
00180   std::cout.setf(ios::scientific);
00181   while ( (!convReached) &&  (nbIter < nbMaxIter) ) {
00182 
00183 #ifdef QC_REPORT_TIME      
00184     timelog << std::endl << "*********** ITER " << nbIter  << " ***********" << std::endl;
00185     timer.start();
00186 #endif 
00187     //
00188     //  Fock est deja construite, on diagonalise
00189     //
00190     fockF.diagonalize(manager, workingSystem, nbIter);
00191     //
00192 #ifdef QC_REPORT_TIME
00193     timer.pause();
00194     diagTimes[nbIter]  = timer.getvalsec();
00195     tDiagTime         += diagTimes[nbIter];
00196     timelog << " * fock diag time:\t"  << diagTimes[nbIter]  << " s" << std::endl;
00197     timer.start();
00198 #endif
00199     //
00200     densityP.computeElems(memory, nbElectrons, eigenVect);
00201     //
00202 #ifdef QC_REPORT_TIME
00203     timer.pause();
00204     densTimes[nbIter] = timer.getvalsec();
00205     tDensTime        += densTimes[nbIter];
00206     timelog << " * density time:\t" << densTimes[nbIter] << " s" << std::endl;
00207     timer.start();
00208 #endif
00209     //
00210     //  Build Fock matrix 
00211     //     Step 1 :  F = H
00212     //
00213     fockF.getMatrix().copy(H);
00214     //     Step 2 :  F = H + 1/2 G(D)
00215     fockF.completeElems(manager, workingSystem, false); 
00216     //
00217 #ifdef QC_REPORT_TIME
00218     timer.pause();
00219     fockTimes[nbIter] = timer.getvalsec();
00220     tFockTime        += fockTimes[nbIter];
00221     timelog << " * fock time:\t\t"  << fockTimes[nbIter] << " s" << std::endl;
00222 #endif
00223     //
00224     //   La nouvelle energie 
00225     //
00226     elecEnergyPrev = elecEnergy;
00227     elecEnergy     = model.computeElecEnergy(workingSystem, densityP);
00228     totalEnergy    = coreEnergy + elecEnergy;
00229     //
00230     //  Test de convergence 
00231     //
00232     error       = D.deltaElemMax(Dprev) ;
00233     convReached = ( (error < threshold) && (QCAbs(elecEnergy - elecEnergyPrev) <= threshold) );
00234 #ifdef QC_VERBOSE_SCF
00235     std::cout << " ** iteration " << nbIter << "/" << nbMaxIter  
00236               << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00237               << std::setprecision(10) << elecEnergy  << std::endl ;
00238 #endif
00239     //
00240     // Si Level Shifting 
00241     //
00242     if(doLevelShifting) {
00243       gap  = eigenVal[nbOccupiedOM] - eigenVal[nbOccupiedOM-1];
00244       gap -= levelShiftingParam;
00245       fockF.levelShifting(workingSystem, gap);
00246     }
00247     Dprev.copy(D); 
00248     ++nbIter;
00249   } 
00250   if(convReached) {
00251     std::cout <<std::endl << "  Convergence in "  ;}
00252   else {
00253     std::cout <<std::endl << "  No convergence in "  ;}
00254   std::cout << nbIter -1 << " iterations (" << nbMaxIter<<")"   
00255             << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00256             << std::setprecision(10) << elecEnergy  << std::endl << std::endl ;
00257 
00258   //
00259   //  La matrice de densite est-elle indempotente ?
00260   //
00261   QCFloat indempo = D.traceProduct(D) - QC_TWO * D.trace() ;
00262   if (convReached  && QCAbs(indempo) > threshold) {
00263     convReached = false;
00264     std::cout << "\n WARNING: The demsity matrix is not indempotent !!\n" << std::endl;
00265   }
00266   //
00267 #ifdef QC_REPORT_TIME
00268   timelog << std::endl << std::endl
00269           << "#CUMULATED TIMES: Pid " << QCCommon::getPidStr() << std::endl << std::endl
00270           << " * fock diag time:\t\t" << tDiagTime << " s" << std::endl
00271           << " * density time:\t\t"   << tDensTime << " s" << std::endl
00272           << " * fock time:\t\t"      << tFockTime << " s" << std::endl;
00273   //
00274   timelog.close();
00275   //
00276   delete [] diagTimes;
00277   delete [] densTimes;
00278   delete [] fockTimes;
00279 #endif
00280 
00281 
00282 #ifndef DEV_SCF_CLASS_NO_MEMORY
00283    Dprev.giveBackTmpMem(memory);
00284 #endif
00285 
00286 }
00287 
00291 template <class TPManager>
00292 void
00293 QCSCFAlgo::optimalDamping (TPManager& manager) {
00294   //
00295 #ifdef QC_VERBOSE_SCF
00296   std::cout << "/******************************************"  << std::endl
00297             << " ********* SOLVE: OPTIMAL DAMPING *********"  << std::endl
00298             << " ******************************************/" << std::endl << std::endl;
00299 #endif
00300 
00301 #ifdef QC_REPORT_TIME
00302   // Les timer de mesures
00303   QCChrono   timer;
00304   QCFloat  * diagTimes       = new QCFloat [nbMaxIter];
00305   QCFloat    tDiagTime       = 0.0;
00306   QCFloat  * densTimes       = new QCFloat [nbMaxIter];
00307   QCFloat    tDensTime       = 0.0;
00308   QCFloat  * fockTimes       = new QCFloat [nbMaxIter];
00309   QCFloat    tFockTime       = 0.0;
00310 
00311   string logname;
00312   logname = ( QCCommon::outdir + PATH_SEPARATOR + "TIME_SCFSolver@" + QCCommon::getPidStr() + OUTPUT_SUFFIX );
00313   ofstream timelog (logname.c_str(), ios::out);  
00314 #endif
00315   //
00316   // le gestionnaire de memoire
00317   //
00318   QCMemory& QCRestrict memory = manager.getMemory();
00319   //
00320   // Le systeme entier
00321   //
00322   typename TPManager::TSystem& QCRestrict system =  manager.getSystem();
00323   //
00324   // Le modele
00325   //
00326   typename TPManager::TModel& model = manager.getModel();
00327   //
00328   // L'iterateur
00329   //
00330   typename TPManager::TSystem::QCIterator workingSystem =  system.begin();
00331   //
00332   // les matrices courantes
00333   //
00334   QCModelMatrices& QCRestrict matrices = workingSystem->getMatrices();
00335   //
00336   // Les alias vers les matrices
00337   //
00338   QCFockGtr<QCSymMatrix>&    QCRestrict fockF     = matrices.getFockF();
00339   QCDensityGtr<QCSymMatrix>& QCRestrict densityP  = matrices.getDensityP();
00340   QCMatElemGtr<QCMatrix>&    QCRestrict eigenVect = matrices.getEigenVectCt();
00341 
00342   QCSymMatrix& QCRestrict D    = densityP.getMatrix();
00343   QCSymMatrix& QCRestrict H    = matrices.getHamiltonH().getMatrix();
00344   QCFloat *QCRestrict eigenVal = matrices.getEigenVal();
00345   //
00346   // Les energies
00347   //
00348   QCFloat coreEnergy, elecEnergy, elecEnergyPrev, totalEnergy;
00349   //
00350   // Les matrices temporaires
00351   //
00352   QCDensityGtr<QCSymMatrix>& QCRestrict Dprev = matrices.getDensityPnm1();
00353   QCDensityGtr<QCSymMatrix>& QCRestrict Dtild = *matrices.getDensityPTild();
00354   QCFockGtr<QCSymMatrix>&    QCRestrict Ftild = *matrices.getFockFTild();
00355 
00356 #ifndef DEV_SCF_CLASS_NO_MEMORY
00357   int dim = D.getDim();
00358   Dprev.sizeMatrix(dim);   Dprev.getMatrix().takeTmpMem(memory); 
00359   Dtild.sizeMatrix(dim);   Dtild.getMatrix().takeTmpMem(memory);
00360   Ftild.sizeMatrix(dim);   Ftild.getMatrix().takeTmpMem(memory);
00361 #endif
00362   //
00363   QCFloat error = 1000000 ;
00364   QCFloat c1,t, s, c, lambdaOpt ;
00365   QCFloat energy_tild, energy_tild_1e, energy_tild_2e, energy_1e ;
00366   //
00367   int nbElectrons   = workingSystem->getNbElectrons();    
00368   int nbOccupiedOM  = nbElectrons / 2;
00369   //
00370   // Construction des matrices
00371   //
00372   system.buildMatrices(manager);
00373   //
00374   // Calcul des energies du systeme  
00375   //
00376   coreEnergy  = model.computeCoreEnergy (manager, workingSystem);
00377   elecEnergy  = model.computeElecEnergy (workingSystem, densityP);
00378   totalEnergy = coreEnergy + elecEnergy; 
00379   workingSystem->setEnergy (totalEnergy);
00380   //
00381   // On commence les Iterations
00382   //
00383   system.prepareIterations(manager);
00384   nbIter      = 0;
00385   convReached = false;
00386   //
00387   t = s = c = lambdaOpt = QC_ZERO;
00388   energy_1e      = QC_HALF * H.traceProduct(D) ;
00389   energy_tild    = elecEnergy;
00390   energy_tild_1e = energy_1e;
00391   //
00392   // On initialise F_tilde avec F
00393   //
00394   Ftild.copyMatrix(fockF);
00395   //
00396   // On initialise P_tilde avec P
00397   //
00398   Dtild.copyMatrix(densityP);
00399   //
00400   //
00401   // **************************************************************************
00402   // ****************** les iterations ****************************************
00403   // **************************************************************************
00404   //
00405   while ( (!convReached) &&  (nbIter < nbMaxIter) ) {
00406     
00407 #ifdef QC_REPORT_TIME      
00408     timelog << std::endl << "*********** ITER " << nbIter << " ***********" << std::endl;
00409     timer.start();
00410 #endif
00411     //
00412     //  Fock est deja construite, on diagonalise
00413     //
00414     Ftild.diagonalize(manager, workingSystem, nbIter);
00415     //
00416 #ifdef QC_REPORT_TIME
00417     timer.pause();
00418     diagTimes[nbIter]  = timer.getvalsec();
00419     tDiagTime         += diagTimes[nbIter];
00420     timelog << " * fock diag time:\t" << diagTimes[nbIter] << " s" << std::endl;
00421     timer.start();
00422 #endif
00423     //    
00424     Dprev.copyMatrix(densityP);
00425     densityP.computeElems(memory, nbElectrons, eigenVect);
00426     //
00427 #ifdef QC_REPORT_TIME
00428     timer.pause();
00429     densTimes[nbIter] = timer.getvalsec();
00430     tDensTime        += densTimes[nbIter];
00431     timelog << " * density time:\t" << densTimes[nbIter] << " s" << std::endl;
00432     timer.start();
00433 #endif
00434     //
00435     //  Build Fock matrix 
00436     //     Step 1 :  F = H
00437     //
00438     fockF.getMatrix().copy(H);
00439     //     Step 2 :  F = H + 1/2 G(D)
00440     fockF.completeElems(manager, workingSystem, false); 
00441     //
00442 #ifdef QC_REPORT_TIME
00443     timer.pause();
00444     fockTimes[nbIter] = timer.getvalsec();
00445     tFockTime        += fockTimes[nbIter];
00446     timelog << " * fock time:\t\t" << fockTimes[nbIter] << " s" << std::endl;
00447 #endif
00448     //
00449     elecEnergyPrev = elecEnergy;
00450     elecEnergy     = model.computeElecEnergy(workingSystem, densityP);
00451     totalEnergy    = coreEnergy + elecEnergy;
00452     energy_1e      = QC_HALF * H.traceProduct(D);
00453     //
00454     //  Test de convergence 
00455     //
00456     error       = densityP.deltaElemMax(Dprev) ;
00457     convReached = ( (error < threshold) && (QCAbs(elecEnergy - elecEnergyPrev) <= threshold) );
00458 #ifdef QC_VERBOSE_SCF
00459     std::cout << " ** iteration " << nbIter << "/" << nbMaxIter  
00460               << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00461               << std::setprecision(10) << elecEnergy  << std::endl ;
00462 #endif
00463     //
00464     // *******************************
00465     // ** Minimisation de l'energie **
00466     // *******************************
00467     //
00468     lambdaOpt = QC_ONE ;
00469     t         = QC_ZERO;
00470     for (int i = 0; i < nbOccupiedOM; ++i) {
00471       t += eigenVal[i];
00472     }
00473     energy_tild_2e  = energy_tild - QC_TWO * energy_tild_1e ;
00474     s               = QC_TWO * t - QC_TWO * (energy_tild_1e + energy_tild_2e);
00475     c               = elecEnergy + energy_tild_2e - QC_TWO * t;
00476     //
00477     // lambdaOpt = min { -s / QCTwo * c), 1}
00478     //
00479     if (c > QC_ZERO && -s < (QC_TWO * c) ) {
00480       lambdaOpt = -s / (QC_TWO * c);
00481     }
00482     c1 = QC_ONE - lambdaOpt ;
00483     //
00484     //   update of Ftild and Dtild
00485     //
00486     Dtild.getMatrix().scale(c1);
00487     Dtild.getMatrix().axpy(lambdaOpt,D);
00488     //
00489     Ftild.getMatrix().scale(c1);
00490     Ftild.getMatrix().axpy(lambdaOpt,fockF.getMatrix());
00491     //
00492     // energy_tild update 
00493     //
00494     energy_tild    += lambdaOpt * ( s + lambdaOpt *c) ;
00495     energy_tild_1e  = c1 * energy_tild_1e + lambdaOpt * energy_1e ;
00496     //
00497     ++nbIter;
00498   }
00499   if(convReached) {
00500     std::cout <<std::endl << "  Convergence in "  ;}
00501   else {
00502     std::cout <<std::endl << "  No convergence in "  ;}
00503   std::cout << nbIter -1 << " iterations (" << nbMaxIter<<")"   
00504             << "   error : "<< std::setprecision(5) <<error<< "  Energie : " 
00505             << std::setprecision(10) << elecEnergy  << std::endl << std::endl ;
00506   
00507   //
00508   //  La matrice de densite est-elle indempotente ?
00509   //
00510   QCFloat indempo = D.traceProduct(D) - QC_TWO * D.trace() ;
00511   if (convReached  && QCAbs(indempo) > threshold) {
00512     convReached = false;
00513     std::cout << "\n WARNING: The demsity matrix is not indempotent !!\n" << std::endl;
00514   }
00515 #ifdef QC_REPORT_TIME
00516   timelog << std::endl << std::endl << "#CUMULATED TIMES: Pid " << QCCommon::getPidStr()
00517           << std::endl << std::endl << " * fock diag time:\t"   << tDiagTime << " s" << std::endl
00518           << " * density time:\t" << tDensTime << " s" << std::endl  << " * fock time:\t\t" 
00519           << tFockTime << " s" << std::endl;
00520   timelog.close();
00521   //
00522   delete [] diagTimes;
00523   delete [] densTimes;
00524   delete [] fockTimes;
00525 #endif
00526   //
00527 #ifndef DEV_SCF_CLASS_NO_MEMORY
00528    Ftild.getMatrix().giveBackTmpMem(memory);
00529    Dtild.getMatrix().giveBackTmpMem(memory);
00530    Dprev.getMatrix().giveBackTmpMem(memory);
00531 #endif
00532 }
00533 
00534 
00535 
00536 
00537 
00538 
00539 QCMANAGER_METH_EXPL_INST_GLOBAL(void QCSCFAlgo::solve);

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