00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00091
00092 QCMemory& QCRestrict memory = manager.getMemory();
00093
00094
00095
00096 TSystem& QCRestrict system = manager.getSystem();
00097
00098
00099
00100
00101 const TParam * QCRestrict parameters = manager.getParameters();
00102
00103
00104
00105 typename TSystem::QCIterator current;
00106
00107
00108
00109 int myrank = system.getRank() ;
00110 bool master = (system.getRank() ==0 );
00111
00112 #ifndef DEV_SCF_CLASS_NO_MEMORY
00113
00114
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
00143
00144 system.buildMatrices(manager);
00145
00146 #ifdef QC_REPORT_TIME
00147
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
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
00178
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
00191
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
00203
00204 system.prepareIterations(manager);
00205 nbIter = 1;
00206 convReached = false;
00207
00208
00209
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
00220
00221 for (current = system.begin(); current != system.end(); current++) {
00222
00223 #ifdef QC_REPORT_TIME
00224 timer.start();
00225 #endif
00226
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
00239 current->getDensityPnm1().swapMatrix(current->getDensityP());
00240 current->computeBFactor(parameters);
00241 }
00242 #ifdef QC_REPORT_TIME
00243 timer.start();
00244 #endif
00245
00246
00247
00248
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
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
00304
00305
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
00318
00319
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
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
00344
00345 if (doLevelShifting) {
00346 std::cerr <<" WARNING LEVEL SHIFTING NOT IMPLENETED IN DIVID AND CONQUER METHOD" <<std::endl;
00347 exit(1) ;
00348
00349
00350 gap -= this->levelShiftingParam;
00351
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
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
00416
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
00438
00439 QCFloat homoOfSD, lumoOfSD;
00440
00441 int nbOccupiedOM, numSD;
00442
00443
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
00452
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
00463
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 }
00473
00474 }
00475 template <class TPManager>
00476 void QCDCAlgo::applyLevelShifting (TPManager& managerInst) {
00477
00478
00479 QCFloat globalGap, globalLumo, globalHomo ;
00480
00481 int numSD;
00482
00483
00484 this->calcHomoLumo(globalLumo, globalHomo);
00485
00486 globalGap = globalLumo - globalHomo;
00487 cout << "globalGap = " << globalGap << endl;
00488
00489
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 }
00496 }
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
00511
00512 QCMemory& QCRestrict memory = manager.getMemory();
00513
00514
00515
00516 TSystem& QCRestrict system = manager.getSystem();
00517
00518
00519
00520 const TParam * QCRestrict parameters = manager.getParameters();
00521
00522
00523
00524 typename TSystem::QCIterator current;
00525
00526
00527
00528 int myrank = system.getRank() ;
00529 bool master = (myrank ==0 );
00530
00531
00532
00533 #ifndef DEV_SCF_CLASS_NO_MEMORY
00534
00535
00536
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
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
00588
00589
00590 QCFloat elecEnergy , elecEnergyPrev ;
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
00598
00599
00600
00601
00602 system.buildMatrices(manager);
00603
00604
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
00623
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
00636
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
00651
00652
00653
00654
00655
00656
00657 nbIter = 0;
00658 convReached = false;
00659 QCFloat error;
00660
00661
00662
00663
00664
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
00674
00675 matrices.getFockFTild()->copyMatrix(matrices.getFockF());
00676 #if defined (QC_TRACE_IN_FILE)
00677 qctrace[ii].writeFockMatrix(matrices);
00678 #endif
00679
00680
00681
00682 matrices.getDensityPTild()->copyMatrix(matrices.getDensityP());
00683 }
00684
00685 nbIter = 1;
00686 convReached = false;
00687
00688
00689
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
00707
00708 for (current = system.begin(); current != system.end(); current++) {
00709
00710 #ifdef QC_REPORT_TIME
00711 timer.start();
00712 #endif
00713
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
00726 current->getDensityPnm1().swapMatrix(current->getDensityP());
00727 current->computeBFactor(parameters);
00728 #if defined (QC_TRACE_IN_FILE)
00729
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
00741
00742
00743 system.adjustFermiEnergySeq(system, memory);
00744
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
00767
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
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
00811
00812
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
00834
00835 elecEnergy = system.completeFockMatricesAndElecEnergy(manager, false, nbIter);
00836
00837
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
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
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
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874 energy_tild += lambdaOpt * ( s + lambdaOpt *c) ;
00875 energy_tild_1e = c1 * energy_tild_1e + lambdaOpt * energy_1e ;
00876 #endif
00877
00878
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
00891
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
00916
00917
00918
00919
00920
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
00944
00945
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