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