00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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 ) {
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
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
00104
00105 QCMemory& QCRestrict memory = manager.getMemory();
00106
00107
00108
00109 typename TPManager::TSystem& QCRestrict system = manager.getSystem();
00110
00111
00112
00113 typename TPManager::TModel& model = manager.getModel();
00114
00115
00116
00117 typename TPManager::TSystem::QCIterator workingSystem = system.begin();
00118
00119
00120
00121 QCModelMatrices& QCRestrict matrices = workingSystem->getMatrices();
00122
00123
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
00134
00135 QCFloat coreEnergy, elecEnergy, elecEnergyPrev, totalEnergy;
00136
00137
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
00154
00155 system.buildMatrices(manager);
00156
00157
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
00171
00172 system.prepareIterations(manager);
00173 nbIter = 0;
00174 convReached = false;
00175
00176
00177
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
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
00211
00212
00213 fockF.getMatrix().copy(H);
00214
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
00225
00226 elecEnergyPrev = elecEnergy;
00227 elecEnergy = model.computeElecEnergy(workingSystem, densityP);
00228 totalEnergy = coreEnergy + elecEnergy;
00229
00230
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
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
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
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
00317
00318 QCMemory& QCRestrict memory = manager.getMemory();
00319
00320
00321
00322 typename TPManager::TSystem& QCRestrict system = manager.getSystem();
00323
00324
00325
00326 typename TPManager::TModel& model = manager.getModel();
00327
00328
00329
00330 typename TPManager::TSystem::QCIterator workingSystem = system.begin();
00331
00332
00333
00334 QCModelMatrices& QCRestrict matrices = workingSystem->getMatrices();
00335
00336
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
00347
00348 QCFloat coreEnergy, elecEnergy, elecEnergyPrev, totalEnergy;
00349
00350
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
00371
00372 system.buildMatrices(manager);
00373
00374
00375
00376 coreEnergy = model.computeCoreEnergy (manager, workingSystem);
00377 elecEnergy = model.computeElecEnergy (workingSystem, densityP);
00378 totalEnergy = coreEnergy + elecEnergy;
00379 workingSystem->setEnergy (totalEnergy);
00380
00381
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
00393
00394 Ftild.copyMatrix(fockF);
00395
00396
00397
00398 Dtild.copyMatrix(densityP);
00399
00400
00401
00402
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
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
00436
00437
00438 fockF.getMatrix().copy(H);
00439
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
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
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
00478
00479 if (c > QC_ZERO && -s < (QC_TWO * c) ) {
00480 lambdaOpt = -s / (QC_TWO * c);
00481 }
00482 c1 = QC_ONE - lambdaOpt ;
00483
00484
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
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
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);