00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef QC_SUBDOMAIN_H_
00022 #define QC_SUBDOMAIN_H_
00023
00024 #include <fstream>
00025 #include <vector>
00026
00027 #include "QCCommon.hpp"
00028 #include "QCGeneralData.hpp"
00029 #include "QCSystem.hpp"
00030 #include "QCModelMatrices.hpp"
00031 #include "QCEnergies.hpp"
00032 #include "QCPartSubDomain.hpp"
00033
00034
00038 #define ATOM_OFFSET 2
00039 #define DOM_ID_IDX 1
00040 #define NB_ATOM_IDX 0
00041
00042
00043
00047 #define OVRLP_MAP_HEAD 1
00048 #define OVRLP_MAP_DECL 1
00049
00050 #define NB_OVRLP_IDX 0
00054 #define IND_MAP_HEAD 7
00055 #define IND_MAP_CELL 3
00059 #define SBD_ID_IDX 0
00060 #define NB_SHARED_ATOM_IDX 1
00061 #define NB_ATOM_SHELL1_IDX 2
00062 #define NB_ATOM_SHELL2_IDX 3
00063 #define NB_SHARED_AO_IDX 4
00064 #define NB_AO_SHELL1_IDX 5
00065 #define NB_AO_SHELL2_IDX 6
00066 //
00067 #define LOCAL_IND_IDX 0
00068 #define REMOTE_IND_IDX 1
00069 #define REMOTE_ZONE_IDX 2
00074
00079 //
00080 class QCSubDomain : public QCSystem, public QCModelMatrices,
00081 public QCEnergies
00082 {
00083
00084 friend ostream&
00085 operator << (ostream& out, const QCSubDomain& domain);
00086
00088 public:
00089
00090
00091
00092
00094 QCSubDomain (void) ;
00098 QCSubDomain (const QCSystem *root);
00099
00100
00101
00102 ~QCSubDomain (void);
00106
00107 int getId (void) const
00108 { return domainId; }
00111 int getNbOverlpSD (void) const
00112 { return overlapMap[NB_OVRLP_IDX]; }
00114 int getGlobalIdx (const int idx) const {
00115 int ret = -1;
00116 if (globalIdx) {
00117 ret = globalIdx[idx];
00118 }
00119 return ret;
00120 }
00121
00122 int getNbAtomCoreShell1 (void) const { return nbAtomCoreShell1; }
00123 int getNbAtomShell2 (void) const { return nbAtomShell2; }
00124 int getNbAOCoreShell1 (void) const { return nbAOCoreShell1; }
00125 int getNbAOShell2 (void) const { return nbAOShell2; }
00126 int getDensityBlockIdx (int idx) const { return densityBlockIdx[idx]; }
00127 int getDensityBlockCount (void) const { return densityBlockCount; }
00129 QCSubDomainZone getZone (int idx) const { return zone[idx]; }
00131 QCSubDomainZone* getZone () const { return zone; }
00132 const int * getGlobalIdxs (void) const { return globalIdx; }
00133 QCFloat * getBFactor (void) { return bFactor; }
00134 QCFloat getBFactor (int idx) const { return bFactor[idx]; }
00135 int getNbFilledOM (void) const { return nbFilledOM; }
00136 int getNbOccupiedOM (void) const { return nbOccupiedOM; }
00137 QCFloat * getOrbitalOccupN (void) { return orbitalOccupN; }
00138 QCFloat getOrbitalOccupN (int idx) const { return orbitalOccupN[idx]; }
00139 int getOverlapMapSize (void) const { return overlapMapSize; }
00140 int * getOverlapMap (void) { return overlapMap; }
00141
00142
00146 void setRootSystem (const QCSystem *root) { rootSystem = root; }
00147 void setZone (int idx, QCSubDomainZone val) { zone[idx] = val; }
00148 void setGlobalIdx (int idx, int val) { globalIdx[idx] = val; }
00149 void setDensityBlockIdx (int idx, int val) { densityBlockIdx[idx] = val; }
00150 void setId (int id) { domainId = id; }
00151 void setGlobalIndex (const int id) { globalDomainId = id; }
00152
00153
00154 void setOrbitalOccupN (int idx, QCFloat val) { orbitalOccupN[idx] = val; }
00155
00156 void setNbFilledOccupiedOM (int filled, int occup) {
00157 nbFilledOM = filled; nbOccupiedOM = occup;
00158 }
00162 QCSystem& getSystem (void) { return QCSystem::asBase(); }
00166 const QCSystem& getSystem (void) const { return QCSystem::asBase(); }
00170 const QCSystem& getRootSystem (void) const { return *rootSystem; }
00174 QCModelMatrices& getMatrices (void) { return QCModelMatrices::asBase(); }
00175
00176
00177 void writeOverlapStructure(){
00178 std::cout <<"*******************************************" <<std::endl;
00179 for(int i =0 ; i <overlapMapSize ; ++i){std::cout << " "<<overlapMap[i] ; }
00180 std::cout <<std::endl <<"*******************************************" <<std::endl;
00181 }
00182
00186 int getNbSharedAtoms (int mapidx) const {
00187 int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00188 return overlapMap[offset + NB_SHARED_ATOM_IDX];
00189 }
00190
00194 int getRemDomainId (int mapidx) const {
00195 int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00196 return overlapMap[offset + SBD_ID_IDX];
00197 }
00201 int getIndLocalIdx (int mapidx, int cellidx) {
00202 int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00203 int cell = offset + IND_MAP_HEAD + cellidx * IND_MAP_CELL;
00204 return overlapMap[cell + LOCAL_IND_IDX];
00205 }
00209 int getIndRemoteIdx (int mapidx, int cellidx) {
00210 int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00211 int cell = offset + IND_MAP_HEAD + cellidx * IND_MAP_CELL;
00212 return overlapMap[cell + REMOTE_IND_IDX];
00213 }
00214
00218 int getMapOffset (int sbdId) const {
00219 int ret = -1;
00220 int offset;
00221 for (int i=0; i < getNbOverlpSD(); ++i) {
00222 offset = overlapMap[OVRLP_MAP_HEAD + i];
00223 if (overlapMap[offset + SBD_ID_IDX] == sbdId) {
00224 ret = offset;
00225 break;
00226 }
00227 }
00228 return ret;
00229 }
00230
00234 int getRemoteIdx (int sbdId, int atomIdx) const {
00235 int idx = -1;
00236 int offset = getMapOffset(sbdId);
00237 if (offset >= 0) {
00238 for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00239 int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00240 if (overlapMap[cellOffset + LOCAL_IND_IDX] == atomIdx) {
00241 idx = overlapMap[cellOffset + REMOTE_IND_IDX];
00242 break;
00243 }
00244 }
00245 }
00246 return idx;
00247 }
00251 bool containsNonShell2 (int mapidx) const {
00252 bool ret = false;
00253 int offset = overlapMap[OVRLP_MAP_HEAD + mapidx];
00254 for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00255 int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00256 if (getZone(overlapMap[cellOffset + LOCAL_IND_IDX]) != QC_SHELL2) {
00257 ret = true;
00258 break;
00259 }
00260 }
00261
00262 return ret;
00263 }
00264
00268 void sortOverlapments (void) {
00269 for (int i=0; i < getNbOverlpSD(); ++i) {
00270 int offset = overlapMap[OVRLP_MAP_HEAD + i];
00271 sort(offset);
00272 }
00273 }
00283 void allocStructures (int nbatoms, int nbao, int nbovrlp, const int * sbdIds, const int * nbshared);
00284
00285 void allocStructures (QCPartSubDomain & domain );
00286
00290 void addOverlpIndirectionInfo (int remoteSbdId,
00291 int atomIdx,
00292 int remoteAtomIdx,
00293 QCSubDomainZone remoteZone);
00297 template <class TPParam>
00298 void setSystemNumbers (const TPParam *QCRestrict parameters);
00302 void fillWeight (const QCGeneralData& data);
00303
00307 template <class TPParam>
00308 void fillDensityArray (const TPParam * QCRestrict parameters,
00309 QCFloat * QCRestrict & densityArray,
00310 bool alloc = false);
00314 template <class TPParam>
00315 void initWeightedDensity (const TPParam * QCRestrict);
00316
00317
00318
00319
00323 template <class TPParam>
00324 void computeBFactor (const TPParam * QCRestrict parameters);
00325
00326
00327
00328
00332 template <class TPParam>
00333 QCFloat traceProductCoreShell1 (const TPParam * QCRestrict parameters,
00334 QCSymMatrix& QCRestrict fockMatrix,
00335 QCSymMatrix& QCRestrict densityMatrix);
00336
00337
00338
00342 void modifyMatricesAndEnergy (const QCFloat& lambdaOpt);
00346 void copyAtomOverlpInto (QCFloat * atomsBuff, int * orderBuff, int *ovrlpBuff,
00347 QCSubDomainZone * zoneBuff, int * blockBuff = NULL) {
00348
00349 atomsBuff[NB_ATOM_IDX] = nbAtoms;
00350 atomsBuff[DOM_ID_IDX] = domainId;
00351 for (int i=0; i < getNbAtoms() * DIMENSION; ++i) {
00352 atomsBuff[i + ATOM_OFFSET] = atoms[i];
00353 }
00354
00355 for (int i=0; i < getNbAtoms() * ORDR_CELL_SIZE; ++i) {
00356 orderBuff[i] = cells[i];
00357 }
00358
00359 for (int i=0; i < getOverlapMapSize(); ++i) {
00360 ovrlpBuff[i] = overlapMap[i];
00361 }
00362
00363 for (int i=0; i < getNbAtoms(); ++i) {
00364 zoneBuff[i] = zone[i];
00365 }
00366
00367 if (blockBuff) {
00368 for (int i=0; i < getNbAtoms(); ++i) {
00369 blockBuff[i] = densityBlockIdx[i];
00370 }
00371 }
00372 }
00376 void setAtomOverlpFrom (QCFloat * atomsBuff, int * orderBuff, int * ovrlpBuff,
00377 QCSubDomainZone * zoneBuff, int * blockBuff = NULL) {
00378
00379 if (atomsBuff) {
00380 nbAtoms = static_cast<int>(atomsBuff[NB_ATOM_IDX]);
00381 domainId = static_cast<int>(atomsBuff[DOM_ID_IDX]);
00382 atoms = atomsBuff + ATOM_OFFSET;
00383 }
00384 cells = orderBuff;
00385 overlapMap = ovrlpBuff;
00386 zone = zoneBuff;
00387 if (blockBuff) {
00388 densityBlockIdx = blockBuff;
00389 }
00390 }
00394 void unsetAtomOverlp (void) {
00395
00396 atoms = NULL;
00397 cells = NULL;
00398 overlapMap = NULL;
00399 zone = NULL;
00400 densityBlockIdx = NULL;
00401 }
00402
00406 void writeInFileNew(std::ofstream& out) ;
00410 void writeSubDomainHeader(std::ofstream & out) ;
00414 void fillAtomInStructure(std::vector<std::vector<QCAtomIn> >& v) ;
00418 void updateInternalStructure() ;
00422 void printHeadOfInternalStructure();
00423
00425 protected:
00426
00427
00429 private:
00430
00431
00435 void swap (int i, int j) {
00436 int tmp;
00437 tmp = overlapMap[i + LOCAL_IND_IDX];
00438 overlapMap[i + LOCAL_IND_IDX] = overlapMap[j + LOCAL_IND_IDX];
00439 overlapMap[j + LOCAL_IND_IDX] = tmp;
00440
00441 tmp = overlapMap[i + REMOTE_IND_IDX];
00442 overlapMap[i + REMOTE_IND_IDX] = overlapMap[j + REMOTE_IND_IDX];
00443 overlapMap[j + REMOTE_IND_IDX] = tmp;
00444
00445 tmp = overlapMap[i + REMOTE_ZONE_IDX];
00446 overlapMap[i + REMOTE_ZONE_IDX] = overlapMap[j + REMOTE_ZONE_IDX];
00447 overlapMap[j + REMOTE_ZONE_IDX] = tmp;
00448 }
00449
00453 void sort (int offset) {
00454 bool sorted = false;
00455 while (!sorted) {
00456 sorted = true;
00457 for (int i=0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]-1; ++i) {
00458 int cell = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00459 if (overlapMap[cell + REMOTE_IND_IDX] >
00460 overlapMap[cell + IND_MAP_CELL + REMOTE_IND_IDX]) {
00461 sorted = false;
00462 swap (cell, cell + IND_MAP_CELL);
00463 }
00464 }
00465 }
00466 }
00467
00470 public:
00471 protected:
00472 private:
00476 int domainId;
00477 int globalDomainId;
00481 const QCSystem * rootSystem;
00485 int nbAtomCoreShell1;
00489 int nbAOCoreShell1;
00493 int nbAtomShell2;
00497 int nbAOShell2;
00502 QCSubDomainZone * zone;
00507 int * globalIdx;
00511 int * densityBlockIdx;
00515 int densityBlockCount;
00519 int * overlapMap;
00523 int overlapMapSize;
00527 QCSymMatrix weight;
00533 QCFloat * bFactor;
00538 int nbFilledOM;
00543 int nbOccupiedOM;
00550 QCFloat * QCRestrict orbitalOccupN;
00551 };
00552
00556 extern ostream& operator << (ostream& out, const QCSubDomain& domain);
00557
00558
00559 #endif // QC_SUBDOMAIN_H_