00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef QC_MD_SYSTEM_H_
00022 #define QC_MD_SYSTEM_H_
00023
00024 #include <vector>
00025
00026 #include "QCSystem.hpp"
00027 #include "QCSubDomain.hpp"
00028 #include "QCEnergies.hpp"
00029 #include "QCPartitioner.hpp"
00030
00031
00032 class QCReader;
00033 template <class TPSystem> class QCPartitioner;
00034
00035
00036
00042
00043
00044
00045 struct QCTopEnergyLevel {
00046
00047
00048
00049
00050 QCTopEnergyLevel (void) {}
00051
00055 void setAttributes (int numSD_, int numOM_) {
00056 this->numSD = numSD_;
00057 this->numOM = numOM_;
00058 }
00059
00060
00064 bool operator < (const QCTopEnergyLevel& other) const {
00065 return (this->energyLevel < other.energyLevel);
00066 }
00067
00068
00071 int numSD;
00072 int numOM;
00073 QCFloat energyLevel;
00074 };
00075
00076
00077
00078
00084 template <class TPSolver>
00085 class QCMDSystem : public QCSystem, public QCEnergies
00086 {
00087
00088 friend class QCReader;
00089
00090
00091
00092
00093
00095 public:
00096
00100 class QCIterator;
00101 friend class QCIterator;
00102
00103 typedef TPSolver TSolver;
00104
00108 QCMDSystem (void);
00109
00113 QCMDSystem (int& , char **& )
00114 { QCMDSystem();}
00115
00119 ~QCMDSystem (void);
00120
00124 int getNbDomains (void) const { return nbDomains; }
00125
00129 QCSubDomain & getSubDomain(const int i) const
00130 { return this->domains[i];}
00134 int getNproc (void) const { return 1; }
00135
00139 int getRank (void) const { return 0; }
00143 QCIterator begin (void) { return QCIterator(this); }
00144
00148 QCIterator end (void) { return QCIterator(this, nbDomains); }
00149
00150
00151
00155 QCIterator getIterator (int idx) { return QCIterator(this, idx); }
00156
00157
00158
00162 QCSystem& getRootSystem (void) { return QCSystem::asBase(); }
00163
00164
00165
00169 const QCSystem& getRootSystem (void) const { return QCSystem::asBase(); }
00170
00174 void addAtom (const QCPoint3D& coords, int type, int globalIdx,
00175 int domainIdx, QCSubDomainZone zone, int atomIdx);
00179 void fillMaps (int nbsbd, const int * domainIdxs, const QCSubDomainZone * zones,
00180 const int * atomIdxs);
00184 void fillMaps ( const std::vector< QCAtomIn > & v ) ;
00193 bool testConvergence (QCFloat & error, QCFloat threshold);
00194
00195 QCFloat computeSumOfEigenValOccupied ();
00196
00200 template <class TPManager>
00201 void init (TPManager& manager, const string& path);
00202
00203
00207 template <class TPSystem>
00208 void allocatePartitioner (QCPartitioner<TPSystem>*& partitioner, const QCGeneralData& data);
00209
00210 template <class TPSystem>
00211 void buildPartitioner (QCPartitioner<TPSystem>*& partitioner,
00212 const QCGeneralData& data, const QCFiles& files,
00213 const int * charge, const int nb);
00214
00215
00216 #ifdef QC_REPORT_TIME
00217
00220 inline void computeAverages (QCFloat& tResTime,QCFloat& tFermiTime, QCFloat& tDensLTime,
00221 QCFloat& tDensCTime, QCFloat& tFockLTime, QCFloat& tFockCTime);
00222 #endif
00223
00227 inline void allocDomains (const int nb, const int offset = 0);
00228
00232 inline void sortSDOverlapments (void);
00233
00237 template <class TPParam>
00238 inline void setDomainsNumbers (const TPParam* QCRestrict params);
00248 inline void allocStructures (int sbdId, int nbatoms, int nbAO,
00249 int nbovlp, int * sbdIds, int * nbshared);
00256 template <class TPManager>
00257 inline void buildMatrices (TPManager& manager);
00261 template <class TPManager>
00262 inline void prepareIterations (TPManager& manager);
00266 template <class TPManager>
00267 void completeHamiltonMatrices (TPManager& manager);
00276 template <class TPManager>
00277 inline void prepareFockMatrices (TPManager& manager, bool isFirstCall);
00290 template <class TPManager>
00291 QCFloat completeFockMatricesAndElecEnergy(TPManager& manager, bool isFirstCall, int iter);
00295 void getGlobalMinMax (QCFloat& , QCFloat& ) {}
00299 void computeElecEnergy() ;
00303 void getTotalSumNbElecs (QCFloat * , int ) {}
00304
00308 QCFloat getTotalSumTrace (void);
00309
00313 void removeLevels (QCFloat& sumNbElecs, QCTopEnergyLevel * topEnergyLevels);
00314
00318 void addLevels (QCFloat& sumNbElecs, const int * nbOM, QCTopEnergyLevel * topEnergyLevels);
00319
00323 template <class TPMDSystem>
00324 void adjustFermiEnergy (TPMDSystem& mdsystem, QCMemory& memory);
00325
00329 template <class TPMDSystem>
00330 void adjustFermiEnergyNew (TPMDSystem& mdsystem, QCMemory& memory);
00331
00334 template <class TPMDSystem>
00335 void adjustFermiEnergySeq (TPMDSystem& mdsystem, QCMemory& memory) {
00336 adjustFermiEnergy(mdsystem, memory);
00337 }
00338
00342 template <class TPManager>
00343 void completeDensityMatrices (TPManager& manager);
00344
00348 template <class TPParam>
00349 QCFloat applyOptimalDamping (const TPParam * QCRestrict params);
00350
00354 template <class TPManager>
00355 void computeEnergyDerivatives (TPManager& ) {}
00356
00366 template <class TPManager>
00367 void writeDensityOnFileAscii(TPManager& manager, const string& path);
00372 template <class TPManager>
00373 void writeLocalDensityOnFilesAscii(TPManager& manager, const string& path);
00377 template <class TPManager>
00378 void writeOrbitalPartition(TPManager& manager, const string& path);
00382 void writeSystemOnFile(const QCFiles & files, const string& path);
00394 void readDensityFromFileAscii(const QCFiles & files, const string& path);
00397 void readSharedDensityFromFileAscii(std::ifstream& data, const std::string &fileName);
00399 void readDistDensityFromFileAscii(std::ifstream& data, const int nproc);
00403 bool checkData(const QCGeneralData& data) ;
00404
00406 protected:
00407
00408
00412 enum QCWay {
00413 ASC,
00414 DESC
00415 };
00416
00417
00424 QCFloat calculateLambdaOptDamp (const QCFloat& globalTraceFP1_PP1OfSD,
00425 const QCFloat& globalTraceF_PP1OfSD,
00426 const QCFloat& globalTraceFP1_POfSD,
00427 const QCFloat& globalTraceF_POfSD);
00428
00432 template <class TPMDSystem>
00433 void computeDichotomyBounds (TPMDSystem& system,
00434 int * nbOM,
00435 QCFloat& bmin,
00436 QCFloat& bmax);
00437
00438
00442 template <class TPMDSystem>
00443 void adjustFermiEnergyDichotomy (TPMDSystem& system,
00444 int * nbOM,
00445 QCTopEnergyLevel * topLevel,
00446 QCFloat sumNE,
00447 QCFloat bmin,
00448 QCFloat bmax,
00449 QCWay way);
00450
00451
00455 template <class TPMDSystem>
00456 void adjustFermiEnergyMonotone (TPMDSystem& system,
00457 int * nbOM,
00458 QCTopEnergyLevel * topLevel);
00459
00460
00461
00462
00463
00465 private:
00466
00467
00470 public:
00471
00472
00477 class QCIterator
00478 {
00480 public:
00481
00485 QCIterator (void) : idx(0), pointer(NULL)
00486 {}
00490 QCIterator (const QCMDSystem * system) : idx(0), pointer(system)
00491 {}
00495 QCIterator (const QCMDSystem * system, int end) : idx(end), pointer(system)
00496 {}
00500 bool operator ++ (void) {
00501 if (pointer->nbDomains > idx) {
00502 idx++;
00503 return true;
00504 }
00505 return false;
00506 }
00510 bool operator -- (void) {
00511 if (idx > 0) {
00512 idx--;
00513 return true;
00514 }
00515 return false;
00516 }
00517
00521 bool operator ++ (int) {
00522 return operator ++ ();
00523 }
00524
00525
00529 bool operator -- (int) {
00530 return operator -- ();
00531 }
00532
00533
00537 QCSubDomain * operator -> (void) {
00538 return &pointer->domains[idx];
00539 }
00540
00541
00545 bool operator == (const QCIterator& iter) const {
00546 return ( (iter.pointer == pointer) && (iter.idx == idx) );
00547 }
00551 bool operator != (const QCIterator& iter) const {
00552 return !(operator == (iter));
00553 }
00554
00556 private:
00557
00558 int idx;
00559
00560 const QCMDSystem * pointer;
00561 };
00562
00563
00564 protected:
00565
00569 int totalNbDomains;
00573 int nbDomains;
00577 QCSubDomain * domains;
00581 int offset;
00585 bool fermiFCall;
00589 QCPartitionType _partitionType ;
00590
00591 private:
00592
00593
00594 };
00598 template <class TPSolver>
00599 inline void
00600 QCMDSystem<TPSolver>::allocDomains (const int nb, const int offset) {
00601
00602 QC_TRACE("BEGIN QCMDSystem<TPSolver>::allocDomains");
00603
00604 if( (this->domains ==NULL ) || (this->domains != NULL && this->nbDomains != nb) )
00605 {
00606 delete [] this->domains ; this->domains = new QCSubDomain [nb];
00607 this->nbDomains = nb;
00608 this->totalNbDomains = nb;
00609 }
00610 for (int i = 0; i < nb; ++i) {
00611 this->domains[i].setId(i + offset);
00612 this->domains[i].setRootSystem(&getRootSystem());
00613 }
00614 QC_TRACE("END QCMDSystem<TPSolver>::allocDomains");
00615 }
00619 template <class TPSolver>
00620 inline void QCMDSystem<TPSolver>::allocStructures (int sbdId, int nbatoms,int nbAO,
00621 int nbovlp, int * sbdIds, int * nbshared) {
00622 QC_TRACE_INIT("BEGIN QCMDSystem<TPSolver>::allocStructures");
00623 domains[sbdId].setGlobalIndex(sbdId);
00624 domains[sbdId].allocStructures(nbatoms, nbAO, nbovlp, sbdIds, nbshared);
00625 QC_TRACE_INIT("END QCMDSystem<TPSolver>::allocStructures");
00626 }
00627
00631 template <class TPSolver>
00632 template <class TPManager>
00633 inline void
00634 QCMDSystem<TPSolver>::buildMatrices (TPManager & manager) {
00635
00636
00637 QCIterator workingSystem = begin();
00638 if(manager.getGeneralData().getInitDensityMatrix() == QC_DENSITY_FROM_FILE ){
00639 this->readDensityFromFileAscii(manager.getFiles(),".");
00640 }
00641 for (workingSystem = this->begin(); workingSystem != this->end(); workingSystem++) {
00642 workingSystem->QCModelMatrices::buildMatrices(manager, workingSystem);
00643 }
00644 }
00648 template <class TPSolver>
00649 template <class TPManager>
00650 inline void
00651 QCMDSystem<TPSolver>::prepareIterations (TPManager & manager) {
00652
00653 QCIterator workingSystem;
00654 for (workingSystem = begin(); workingSystem != end(); workingSystem++) {
00655 workingSystem->QCModelMatrices::prepareIterations(manager);
00656 }
00657 }
00658
00662 template <class TPSolver>
00663 template <class TPManager>
00664 inline void
00665 QCMDSystem<TPSolver>::prepareFockMatrices (TPManager& manager, bool isFirstCall) {
00666
00667 QCIterator iter;
00668 for (iter = begin(); iter != end(); iter++) {
00669 iter->getFockF().copyMatrix(iter->getHamiltonH());
00670 iter->getFockF().completeElems(manager, iter, isFirstCall);
00671 }
00672 }
00673
00674
00675
00676
00677 #ifdef QC_REPORT_TIME
00678
00681 template <class TPSolver>
00682 inline void
00683 QCMDSystem<TPSolver>::computeAverages (QCFloat& tResTime,
00684 QCFloat& ,
00685 QCFloat& tDensLTime,
00686 QCFloat& ,
00687 QCFloat& tFockLTime,
00688 QCFloat& )
00689 {
00690 tResTime /= nbDomains;
00691 tDensLTime /= nbDomains;
00692 tFockLTime /= nbDomains;
00693 }
00694 #endif
00695
00696
00700 template <class TPSolver>
00701 extern ostream&
00702 operator << (ostream& out, const QCMDSystem<TPSolver>& system);
00703
00704
00705
00706
00707
00708
00709 #endif // QC_MD_SYSTEM_H_