00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <iomanip>
00021 #include <vector>
00022
00023 #include "QCSubDomain.hpp"
00024 #include "QCTable.hpp"
00025 #include "QCMndoParam.hpp"
00026 #include "QCAm1Param.hpp"
00027 #include "QCPm3Param.hpp"
00028
00029
00030
00034 QCSubDomain::QCSubDomain (const QCSystem *root) : QCSystem(), QCModelMatrices(),
00035 QCEnergies(), domainId(-1), globalDomainId(-1),rootSystem(root),
00036 nbAtomCoreShell1(0), nbAOCoreShell1(0), nbAtomShell2(0),
00037 nbAOShell2(0), zone(NULL), globalIdx(NULL), densityBlockIdx(NULL),
00038 densityBlockCount(0), overlapMap(NULL), overlapMapSize(0), weight(),
00039 bFactor(NULL), nbFilledOM(0), nbOccupiedOM(0), orbitalOccupN(NULL)
00040 {}
00041
00042
00043
00044
00048 QCSubDomain::QCSubDomain (void) : QCSystem(), QCModelMatrices(), QCEnergies(), domainId(-1),globalDomainId(-1),
00049 rootSystem(NULL), nbAtomCoreShell1(0), nbAOCoreShell1(0),
00050 nbAtomShell2(0), nbAOShell2(0), zone(NULL), globalIdx(NULL),
00051 densityBlockIdx(NULL), densityBlockCount(0), overlapMap(NULL),
00052 overlapMapSize(0), weight(), bFactor(NULL), nbFilledOM(0),
00053 nbOccupiedOM(0), orbitalOccupN(NULL)
00054 {}
00055
00056
00057
00058
00062 QCSubDomain::~QCSubDomain (void) {
00063 QC_TRACE_END("BEGIN QCSubDomain::~QCSubDomain");
00064 if (zone) {
00065 delete [] zone; zone = NULL;
00066 }
00067 if (globalIdx) {
00068 delete [] globalIdx; globalIdx = NULL;
00069 }
00070 if (densityBlockIdx) {
00071 delete [] densityBlockIdx; densityBlockIdx = NULL;
00072 }
00073 if (overlapMap) {
00074 delete [] overlapMap; overlapMap = NULL;
00075 }
00076 if (bFactor) {
00077 delete [] bFactor; bFactor = NULL;
00078 }
00079 if (orbitalOccupN) {
00080 delete [] orbitalOccupN; orbitalOccupN = NULL;
00081 }
00082 QC_TRACE_END("END QCSubDomain::~QCSubDomain");
00083 }
00084
00085
00086
00087
00091 void
00092 QCSubDomain::allocStructures (int nbatoms, int nbao, int nbovrlp,
00093 const int *sbdIds, const int *nbshared) {
00094 QC_TRACE_INIT("BEGIN QCSubDomain::allocStructures for subdomain "<< domainId);
00095 QC_TRACE_INIT(" nbatoms "<<nbatoms <<" nbao "<< nbao<<" nbovrlp "<<nbovrlp);
00096
00097 allocateAtoms(nbatoms);
00098 allocateCells(nbatoms);
00099 zone = new QCSubDomainZone [nbatoms];
00100 globalIdx = new int [nbatoms];
00101 densityBlockIdx = new int [nbatoms];
00102 bFactor = new QCFloat [nbao];
00103 orbitalOccupN = new QCFloat [nbao];
00104 weight.setDimAndAllocate(nbatoms);
00105
00106 this->overlapMapSize = OVRLP_MAP_HEAD + nbovrlp * OVRLP_MAP_DECL;
00107
00108 int * offset = new int [nbovrlp];
00109
00110 for (int i=0; i < nbovrlp; ++i) {
00111 offset[i] = this->overlapMapSize;
00112 this->overlapMapSize += IND_MAP_HEAD + IND_MAP_CELL * nbshared[i];
00113 }
00114
00115 this->overlapMap = new int [this->overlapMapSize];
00116 this->overlapMap[NB_OVRLP_IDX] = nbovrlp;
00117
00118 for (int i=0; i < nbovrlp; ++i) {
00119
00120 this->overlapMap [OVRLP_MAP_HEAD + i ] = offset[i];
00121 this->overlapMap [offset[i] + SBD_ID_IDX ] = sbdIds[i];
00122 this->overlapMap [offset[i] + NB_SHARED_ATOM_IDX] = 0;
00123 this->overlapMap [offset[i] + NB_ATOM_SHELL1_IDX] = 0;
00124 this->overlapMap [offset[i] + NB_ATOM_SHELL2_IDX] = 0;
00125 this->overlapMap [offset[i] + NB_SHARED_AO_IDX ] = 0;
00126 this->overlapMap [offset[i] + NB_AO_SHELL1_IDX ] = 0;
00127 this->overlapMap [offset[i] + NB_AO_SHELL2_IDX ] = 0;
00128 }
00129
00130 delete [] offset;
00131 QC_TRACE_INIT("END QCSubDomain::allocStructures");
00132 }
00133
00139 void QCSubDomain::allocStructures (QCPartSubDomain & domain ){
00140 QC_TRACE_INIT("BEGIN QCSubDomain::allocStructures (QCPartSubDomain & domain ) ");
00141
00142 int nbAtoms = domain.getNbAtoms(), nbNeighbors= domain.getNbNeighbors(), nbAo = domain.getLoad();
00143
00144 allocateAtoms(nbAtoms);
00145 allocateCells(nbAtoms);
00146
00147 zone = new QCSubDomainZone [nbAtoms];
00148 globalIdx = new int [nbAtoms];
00149 densityBlockIdx = new int [nbAtoms];
00150
00151 weight.setDimAndAllocate(nbAtoms);
00152
00153 bFactor = new QCFloat [nbAo];
00154 orbitalOccupN = new QCFloat [nbAo];
00155
00156 this->overlapMapSize = OVRLP_MAP_HEAD + nbNeighbors * OVRLP_MAP_DECL;
00157 int * offset = new int [nbNeighbors];
00158
00159 for(int i = 0 ; i < domain.getNbNeighbors(); ++i){
00160 offset[i] = this->overlapMapSize;
00161 this->overlapMapSize += IND_MAP_HEAD + IND_MAP_CELL * ( domain.getNbAtomsShell1(i) + domain.getNbAtomsShell2(i));
00162 }
00163 this->overlapMap = new int [this->overlapMapSize];
00164 this->overlapMap[NB_OVRLP_IDX] = nbNeighbors;
00165
00166 for (int i=0; i < nbNeighbors; ++i) {
00167 this->overlapMap [OVRLP_MAP_HEAD + i ] = offset[i];
00168 this->overlapMap [offset[i] + SBD_ID_IDX ] = domain.getNeighborIdx(i) ;
00169 this->overlapMap [offset[i] + NB_SHARED_ATOM_IDX] = 0 ;
00170 this->overlapMap [offset[i] + NB_ATOM_SHELL1_IDX] = domain.getNbAtomsShell1(i);
00171 this->overlapMap [offset[i] + NB_ATOM_SHELL2_IDX] = domain.getNbAtomsShell2(i);
00172 this->overlapMap [offset[i] + NB_SHARED_AO_IDX ] = 0;
00173 this->overlapMap [offset[i] + NB_AO_SHELL1_IDX ] = 0;
00174 this->overlapMap [offset[i] + NB_AO_SHELL2_IDX ] = 0;
00175 }
00176
00177 delete [] offset;
00178
00179 #ifdef DEBUG_INIT_SUBDOMAIN
00180 this->printHeadOfInternalStructure();
00181 #endif
00182 QC_TRACE_INIT("END QCSubDomain::allocStructures (QCPartSubDomain & domain ) ");
00183 }
00184
00188 void
00189 QCSubDomain::addOverlpIndirectionInfo (int remoteSbdId, int atomIdx,
00190 int remoteAtomIdx, QCSubDomainZone remoteZone) {
00191
00192
00193 int offset = this->getMapOffset(remoteSbdId);
00194
00195 #ifndef QC_NO_DEBUG
00196 assert(offset >= 0);
00197 #endif
00198
00199
00200 int& nbCells = overlapMap[offset + NB_SHARED_ATOM_IDX];
00201 int cellOffset = ( offset + IND_MAP_HEAD + nbCells * IND_MAP_CELL);
00202
00203 overlapMap [cellOffset + LOCAL_IND_IDX ] = atomIdx;
00204 overlapMap [cellOffset + REMOTE_IND_IDX ] = remoteAtomIdx;
00205 overlapMap [cellOffset + REMOTE_ZONE_IDX] = remoteZone;
00206
00207 nbCells++;
00208 }
00209
00213 template <class TPParam>
00214 void
00215 QCSubDomain::setSystemNumbers (const TPParam *QCRestrict parameters) {
00216 QC_TRACE_INIT("BEGIN QCSubDomain::setSystemNumbers ");
00217
00218 nbAtomicOrbitals = densityBlockCount = nbElectrons = 0 ;
00219 nbAtomShell2 = nbAtomCoreShell1 = 0 ;
00220
00221 for (int i=0; i < nbAtoms; i++) {
00222
00227 setFirstAOAt(i, nbAtomicOrbitals);
00228 setDensityBlockIdx(i, densityBlockCount);
00229 int nbAO = parameters[getParamIndexAt(i)].getNbAO();
00230 nbAtomicOrbitals += nbAO;
00231 densityBlockCount += (nbAO * (nbAO + 1)) / 2;
00232 nbElectrons += parameters[getParamIndexAt(i)].getNbValenceElec();
00233
00234
00235
00236
00237
00238
00239
00240
00241 }
00245 nbElectrons -= systemCharge * (nbElectrons / rootSystem->getNbElectrons());
00246 spinMultiplicity = rootSystem->getSpinMultiplicity();
00247 #if defined TRACE_OUT
00248 this->printHeadOfInternalStructure();
00249 #endif
00250 QC_TRACE_INIT("END QCSubDomain::setSystemNumbers ");
00251 }
00252
00256 void QCSubDomain::fillWeight (const QCGeneralData& data) {
00257
00258 int k, i, j;
00259 int localIdxA, localIdxB;
00260
00261 for (i = 0; i < getNbAtoms(); ++i) {
00262 for (j = 0; j <= i; ++j) {
00263 weight[i][j] = 1;
00264 }
00265 }
00266
00267
00268 for (k = 0; k < getNbOverlpSD(); ++k) {
00269 int offset = overlapMap[OVRLP_MAP_HEAD + k];
00270
00271
00272 for (i = 0; i < overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00273 int cellA = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00274 localIdxA = overlapMap[cellA + LOCAL_IND_IDX];
00275
00276
00277
00278 for (j = 0; j <= i; ++j) {
00279 int cellB = offset + IND_MAP_HEAD + j * IND_MAP_CELL;
00280 localIdxB = overlapMap[cellB + LOCAL_IND_IDX];
00281
00282
00283 if (localIdxA > localIdxB) {
00284 weight[localIdxA][localIdxB] += 1;
00285
00286 } else {
00287 weight[localIdxB][localIdxA] += 1;
00288 }
00289 }
00290 }
00291 }
00292
00293
00294 for (i = 0; i < getNbAtoms(); ++i) {
00295 for (j = 0; j <= i; ++j) {
00296
00297
00298 if (data.getPartitionType() == QC_DIXON_PART) {
00299
00300 if (zone[i] == QC_SHELL2 || zone[j] == QC_SHELL2) {
00301 weight[i][j] = QC_ZERO;
00302
00303 } else if (zone[i] == QC_SHELL1 && zone[j] == QC_SHELL1) {
00304 weight[i][j] = QC_ZERO;
00305
00306 } else {
00307 weight[i][j] = QC_ONE / weight[i][j];
00308 }
00309
00310 } else if (data.getPartitionType() == QC_STANDARD_PART) {
00311
00312 if (zone[i] == QC_SHELL2 || zone[j] == QC_SHELL2) {
00313 weight[i][j] = QC_ZERO;
00314
00315 } else {
00316 weight[i][j] = QC_ONE / weight[i][j];
00317 }
00318
00319 }
00320 }
00321 }
00322 }
00323
00327 template <class TPParam>
00328 void
00329 QCSubDomain::fillDensityArray (const TPParam * QCRestrict parameters,
00330 QCFloat * QCRestrict & densityArray,
00331 bool alloc) {
00332
00333
00334 const TPParam * QCRestrict param;
00335
00336
00337 const QCSymMatrix& density = getDensityP().getMatrix();
00338
00339
00340 int i, mu, nu, firstAO, lastAO, blockIdx;
00341
00342 if (alloc) {
00343 densityArray = new QCFloat [densityBlockCount];
00344 }
00345
00346 blockIdx = 0;
00347
00348 for (i=0; i < nbAtoms; ++i) {
00349 param = ¶meters[getParamIndexAt(i)];
00350 firstAO = getFirstAOAt(i);
00351 lastAO = firstAO + param->getNbAO();
00352
00353 for (mu = firstAO; mu < lastAO; ++mu) {
00354 for (nu = firstAO; nu <= mu; ++nu) {
00355 densityArray [blockIdx] = density [mu][nu];
00356 ++blockIdx;
00357 }
00358 }
00359 }
00360 }
00361
00367 template <class TPParam>
00368 void
00369 QCSubDomain::initWeightedDensity (const TPParam * QCRestrict parameters) {
00370
00371
00372 QCSymMatrix& QCRestrict densityP = getDensityP().getMatrix();
00373
00374
00375 QCSymMatrix& QCRestrict weightedDensityP =
00376 getWeightedDensityP()->getMatrix();
00377
00378
00379 const TPParam * QCRestrict paramA;
00380 const TPParam * QCRestrict paramB;
00381
00382
00383 int firstAOofA, firstAOofB;
00384 int lastAOofA, lastAOofB;
00385
00386
00387 QCFloat weightAA, weightAB;
00388
00389
00390 int mu, nu, lda;
00391 int a, b;
00392
00393
00394 for (a = 0; a < getNbAtoms(); ++a) {
00395
00396 paramA = ¶meters[getParamIndexAt(a)];
00397 firstAOofA = getFirstAOAt(a);
00398 lastAOofA = firstAOofA + paramA->getNbAO();
00399
00400 for (b = 0; b < a; ++b) {
00401
00402 paramB = ¶meters[getParamIndexAt(b)];
00403 firstAOofB = getFirstAOAt(b);
00404 lastAOofB = firstAOofB + paramB->getNbAO();
00405
00406 weightAB = weight[a][b];
00407
00408 for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00409 for (lda = firstAOofB; lda < lastAOofB; ++lda) {
00410 weightedDensityP[mu][lda] = weightAB * densityP[mu][lda];
00411 }
00412 }
00413 }
00414
00415
00416 weightAA = weight[a][a];
00417
00418 for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00419 for (nu = firstAOofA; nu <= mu; ++nu) {
00420 weightedDensityP[mu][nu] = weightAA * densityP[mu][nu];
00421 }
00422 }
00423 }
00424 }
00425
00429 template <class TPParam>
00430 void
00431 QCSubDomain::computeBFactor (const TPParam * QCRestrict parameters) {
00432
00433
00434 QCMatrix& QCRestrict localEigenVectCt = getEigenVectCt().getMatrix();
00435
00436
00437 int firsAOofB, lastAOofB;
00438
00439
00440 const TPParam * QCRestrict paramB;
00441
00442
00443 QCFloat weightBB;
00444
00445
00446 int b, i, mu;
00447
00448
00449 for (i = 0; i < getNbAtomicOrbitals(); ++i) {
00450 bFactor[i] = QC_ZERO;
00451 for (b = 0; b < getNbAtoms(); ++b) {
00452 weightBB = weight[b][b];
00453 paramB = ¶meters[getParamIndexAt(b)];
00454 firsAOofB = getFirstAOAt(b);
00455 lastAOofB = firsAOofB + paramB->getNbAO();
00456
00457 for (mu = firsAOofB; mu < lastAOofB; ++mu) {
00466 bFactor[i] += weightBB * FQCPow<2>(localEigenVectCt[i][mu]);
00467 }
00468 }
00469
00470 }
00471 }
00472
00473
00474
00475
00476
00477
00478
00479 template <class TPParam>
00480 QCFloat
00481 QCSubDomain::traceProductCoreShell1 (const TPParam * QCRestrict parameters,
00482 QCSymMatrix& QCRestrict fockMatrix,
00483 QCSymMatrix& QCRestrict densityMatrix) {
00484
00485
00486
00487 const TPParam * QCRestrict parameterA;
00488 const TPParam * QCRestrict parameterB;
00489
00490
00491
00492
00493
00494 int firstAOofA, firstAOofB, lastAOofA, lastAOofB;
00495
00496
00497 register QCFloat interVal, interVal2, interVal3;
00498
00499
00500 int i, j, mu, nu, lda;
00501
00502
00503 QCFloat traceFP = QC_ZERO;
00504
00505 for (i = 0; i < getNbAtoms(); ++i) {
00506
00507 if (getZone(i) != QC_SHELL2) {
00508 parameterA = ¶meters[getParamIndexAt(i)];
00509 firstAOofA = getFirstAOAt(i);
00510 lastAOofA = firstAOofA + parameterA->getNbAO();
00511
00512 for (j = 0; j < i; ++j) {
00513
00514 if (getZone(j) == QC_SHELL2) {
00515 continue ;
00516 }
00517
00518 parameterB = ¶meters[getParamIndexAt(j)];
00519 firstAOofB = getFirstAOAt(j);
00520 lastAOofB = firstAOofB + parameterB->getNbAO();
00521
00522
00523
00524 interVal = QC_ZERO;
00525 for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00526 for (lda = firstAOofB; lda < lastAOofB; ++lda) {
00527
00528
00529
00530
00531
00532 interVal += (fockMatrix[mu][lda] * densityMatrix[mu][lda]);
00533 }
00534 }
00535
00536
00537
00538 if (getZone(i) == QC_CORE || getZone(j) == QC_CORE) {
00539 traceFP += interVal;
00540 }
00541
00542
00543 else {
00544
00545 traceFP += weight[i][j] * interVal;
00546 }
00547 }
00548
00549
00550 interVal = QC_ZERO;
00551 interVal2 = QC_ZERO;
00552 for (mu = firstAOofA; mu < lastAOofA; ++mu) {
00553 for (nu = firstAOofA; nu < mu; ++nu) {
00554 interVal += (fockMatrix[mu][nu] * densityMatrix[mu][nu]);
00555 }
00556 interVal2 += (fockMatrix[mu][mu] * densityMatrix[mu][mu]);
00557 }
00558
00559
00560
00561 interVal3 = interVal + QC_HALF * interVal2;
00562
00563
00564
00565 if (getZone(i) == QC_CORE) {
00566 traceFP += interVal3;
00567 }
00568
00569
00570 else {
00571
00572 traceFP += weight[i][i] * interVal3;
00573 }
00574 }
00575 }
00576 return traceFP;
00577 }
00578
00582 void QCSubDomain::modifyMatricesAndEnergy (const QCFloat& lambdaOpt) {
00583
00584 QCFloat energyFockTild;
00585 QCFloat energyTild1e, c1;
00586
00587 c1 = QC_ONE - lambdaOpt ;
00588 fockFTild->getMatrix().scale(c1);
00589 fockFTild->getMatrix().axpy(lambdaOpt, fockF.getMatrix());
00590
00591 densityPTild->getMatrix().scale(c1);
00592 densityPTild->getMatrix().axpy(lambdaOpt, densityP.getMatrix());
00593
00594
00595 fockF.copyMatrix(*fockFTild);
00596
00597
00598
00599 energyFockTild = ( QC_HALF *
00600 fockFTild->getMatrix().traceProduct(densityPTild->getMatrix()) );
00601 energyTild1e = ( QC_HALF *
00602 hamiltonH.getMatrix().traceProduct(densityPTild->getMatrix()) );
00603
00604 elecEnergyETild = energyTild1e + energyFockTild;
00605 elecEnergyETild1e = energyTild1e;
00606 }
00607
00611 ostream& operator << (ostream& out, const QCSubDomain& domain) {
00612
00613 out << " * Atoms:" << endl;
00614 out << "{" << endl << " LOCAL/GLOBAL:" << " ( X"
00615 << ", Y, Z), TYPE, ZONE, (SD i: LOCAL)"
00616 << endl << endl;
00617
00618 for (int i=0; i < domain.getNbAtoms(); ++i) {
00619 QCPoint3D current;
00620
00621 domain.getPointAt(i, current);
00622
00623 out << " " << setw(5) << i << "/ "
00624 << setw(5) << domain.getGlobalIdx(i) << ": "
00625 << current << ", "
00626 << table[domain.getType(i)] << ", "
00627 << QCCommon::domainZones.find(domain.getZone(i))->second
00628 << ", ";
00629
00630 for (int j=0; j < domain.overlapMap[NB_OVRLP_IDX]; ++j) {
00631 int offset = domain.overlapMap[OVRLP_MAP_HEAD + j];
00632 int domainId = domain.overlapMap[offset + SBD_ID_IDX];
00633 int remoteIdx = domain.getRemoteIdx(domainId, i);
00634
00635 if (remoteIdx >= 0) {
00636 out << "(sd" << setw(2) << domainId << ": " << setw(5) << remoteIdx << ") ";
00637 }
00638 }
00639
00640 out << endl;
00641 }
00642 out << "}" << endl;
00643 out << " * Overlapment:" << endl;
00644 out << "{" << endl;
00645
00646 for (int j=0; j < domain.overlapMap[NB_OVRLP_IDX]; ++j) {
00647 int offset = domain.overlapMap[OVRLP_MAP_HEAD + j];
00648
00649 out << "* domain " << setw(2) << domain.getId()
00650 << " <-> domain " << setw(2) << domain.overlapMap[offset + SBD_ID_IDX]
00651 << " (" << setw(4) << domain.overlapMap[offset + NB_SHARED_ATOM_IDX]
00652 << " atoms):" << endl
00653 << " {" << endl;
00654
00655 out << " (LOCAL, GlOBAL ZONE) -> (REMOTE, RZONE)" << endl << endl;
00656 for (int i=0; i < domain.overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00657
00658 int cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00659 int localIdx = domain.overlapMap[cellOffset + LOCAL_IND_IDX];
00660 int remoteIdx = domain.overlapMap[cellOffset + REMOTE_IND_IDX];
00661 int remoteZone = domain.overlapMap[cellOffset + REMOTE_ZONE_IDX];
00662
00663 out << " (" << setw(5) << localIdx << " "<< domain.globalIdx[localIdx]
00664 << ", " << QCCommon::domainZones.find(domain.getZone(localIdx))->second
00665 << ") -> ( " << setw(5) << remoteIdx << ", "
00666 << QCCommon::domainZones.find(static_cast<QCSubDomainZone>(remoteZone))->second
00667 << ")";
00668 out << endl;
00669 }
00670 out << " }" << endl << endl;
00671 }
00672 out << "}" << endl;
00673
00674 out << "\n * Nb Atoms: " << domain.getNbAtoms()
00675 << " (core+shell1: " << domain.getNbAtomCoreShell1()
00676 << ", shell2: " << domain.getNbAtomShell2()
00677 << ")"
00678 << "\n * System Charge: " << domain.getSystemCharge()
00679 << "\n * Spin Multiplicity: " << domain.getSpinMultiplicity()
00680 << "\n * NbAtomicOrbitals: " << domain.getNbAtomicOrbitals()
00681 << " (core+shell1: " << domain.getNbAOCoreShell1()
00682 << ", shell2: " << domain.getNbAOShell2()
00683 << ")"
00684 << "\n * NbElectrons: " << domain.getNbElectrons()
00685 << endl;
00686
00687 return out;
00688 }
00689
00690 void QCSubDomain::writeInFileNew(std::ofstream& out){
00691
00692 QC_TRACE_OUT("BEGIN QCSubDomain::writeInFileNew");
00693
00694 #if defined TRACE_OUT
00695 this->printHeadOfInternalStructure();
00696 #endif
00697 int nbSD, numCore, numOverlap, numCut, n, j;
00698
00699 nbSD = this->overlapMap[NB_OVRLP_IDX];
00700 numCore = numOverlap = numCut = 0 ;
00702 for ( n = 0 ; n < this->getNbAtoms() ; ++n){
00703 if (this->zone[n] == QC_CORE){
00704 ++numCore ; }
00705 else if (this->zone[n] == QC_SHELL1){
00706 ++numOverlap ; }
00707 else {
00708 ++numCut ; }
00709 }
00710 out << "Sub Domain number : " << this->domainId <<std::endl ;
00711 out << " "<< nbSD <<" " << numCore <<" "<< numOverlap << " " << numCut<<std::endl ;
00712 out << " CORE ZONE : " <<std::endl ;
00713 if (numCore !=0){
00714 out << " ";
00715 for ( n = 0 ; n < this->getNbAtoms() ; ++n){
00716 if (this->zone[n] == QC_CORE){
00717 out << this->globalIdx[n] << " " ;
00718 }
00719 }
00720 out << std::endl ;
00721 }
00722 #ifdef DEBUG_INIT_SUBDOMAIN
00723 this->printHeadOfInternalStructure();
00724 #endif
00725 int offset , nbShared, nAt ,localIdx;
00726 out << " OVERLAP ZONE : " <<std::endl ;
00727 for ( n = 0 ; n < nbSD ; ++n){
00728 offset = this->overlapMap[OVRLP_MAP_HEAD + n*OVRLP_MAP_DECL] ;
00729 nbShared = this->overlapMap[offset + NB_SHARED_ATOM_IDX] ;
00730 nAt = this->overlapMap[offset + NB_ATOM_SHELL1_IDX] ;
00731 if (nAt == 0 ) { continue;}
00732 out << " " << this->overlapMap[offset + SBD_ID_IDX] <<" "<< nAt << " ";
00733
00734 offset += IND_MAP_HEAD + LOCAL_IND_IDX;
00735 for( j = 0 ; j < nbShared ; ++j){
00736 localIdx = overlapMap[offset + j*IND_MAP_CELL];
00737 if( this->zone[localIdx] == QC_SHELL1){
00738 out << this->globalIdx[localIdx]<< " " ;
00739 }
00740 }
00741 out << std::endl;
00742 }
00743
00744 out << " CUT ZONE : " <<std::endl ;
00745 for ( n = 0 ; n < nbSD ; ++n){
00746 offset = this->overlapMap[OVRLP_MAP_HEAD + n*OVRLP_MAP_DECL] ;
00747 nbShared = this->overlapMap[offset + NB_SHARED_ATOM_IDX] ;
00748 nAt = this->overlapMap[offset + NB_ATOM_SHELL2_IDX] ;
00749
00750 if (nAt == 0 ) { continue;}
00751
00752 out << " " << this->overlapMap[offset + SBD_ID_IDX] <<" "<< nAt << " ";
00753 offset += IND_MAP_HEAD + LOCAL_IND_IDX ;
00754 for( j = 0; j < nbShared ; ++j){
00755 localIdx = overlapMap[offset + j*IND_MAP_CELL];
00756
00757
00758 if( this->zone[localIdx] == QC_SHELL2){
00759 out << this->globalIdx[localIdx]<< " " ;
00760 }
00761 }
00762 out << std::endl;
00763 }
00764
00765 QC_TRACE_OUT("END QCSubDomain::writeInFileNew");
00766 }
00767
00769 void QCSubDomain::updateInternalStructure(){
00771 int n,i, offset, localIdx,cellOffset, nbSHELL1, nbSHELL2, nbOASHELL1, nbOASHELL2;
00773 QC_TRACE_INIT("BEGIN QCSubDomain::updateInternalStructure "<< this->domainId << " GlobalIdx " <<globalDomainId);
00774 #ifdef DEBUG_INIT_SUBDOMAIN
00775 std::cout<< std::endl<< std::endl<< " HEAD updated "<< std::endl
00776 <<" Number Of SDs : "<<this->overlapMap[NB_OVRLP_IDX]<< std::endl;
00777 this->printHeadOfInternalStructure();
00778 #endif
00779
00780
00781 nbOASHELL1 = nbOASHELL2 = 0 ;
00782 int oa;
00784
00785 for ( n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00786
00787 offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00789 nbSHELL1 = nbSHELL2 = 0 ; nbOASHELL1 = nbOASHELL2 = 0 ;
00790 for ( i=0; i < this->overlapMap[offset + NB_SHARED_ATOM_IDX]; ++i) {
00791
00792 cellOffset = offset + IND_MAP_HEAD + i * IND_MAP_CELL;
00793 localIdx = this->overlapMap[cellOffset + LOCAL_IND_IDX];
00794
00795 oa = (this->QCSystem::getType(localIdx) == 1 )? 1 : 4 ;
00796 if (this->zone[localIdx] == QC_SHELL1){
00797 ++nbSHELL1 ; nbOASHELL1 += oa ;
00798 }
00799 else if (this->zone[localIdx] == QC_SHELL2){
00800 ++nbSHELL2; nbOASHELL2 += oa ;
00801 }
00802 }
00803 this->overlapMap[offset + NB_ATOM_SHELL1_IDX] = nbSHELL1 ;
00804 this->overlapMap[offset + NB_AO_SHELL1_IDX] = nbOASHELL1 ;
00805 this->overlapMap[offset + NB_ATOM_SHELL2_IDX] = nbSHELL2 ;
00806 this->overlapMap[offset + NB_AO_SHELL2_IDX] = nbOASHELL2 ;
00807 this->overlapMap[offset + NB_SHARED_ATOM_IDX] = nbSHELL1+ nbSHELL2;
00808 this->overlapMap[offset + NB_SHARED_AO_IDX] = nbOASHELL1+ nbOASHELL2;
00809 }
00810
00811
00812
00813 this->nbAtomCoreShell1 = this->nbAOCoreShell1 = this->nbAtomShell2 = this->nbAOShell2 = 0 ;
00814 for (int at =0; at < this->QCSystem::getNbAtoms() ;++at){
00815
00816 oa = (this->QCSystem::getType(at) == 1 )? 1 : 4 ;
00817 if (this->zone[at] == QC_SHELL2){
00818 ++this->nbAtomShell2 ; this->nbAOShell2 += oa;
00819 }
00820 else {
00821 this->nbAOCoreShell1 += oa; ++this->nbAtomCoreShell1;
00822 }
00823 }
00824 #ifdef DEBUG_INIT_SUBDOMAIN
00825 this->printHeadOfInternalStructure();
00826 #endif
00827
00828 QC_TRACE_INIT("END QCSubDomain::updateInternalStructure");
00829 }
00830 void QCSubDomain::printHeadOfInternalStructure(){
00832 int n, offset;
00836 QC_TRACE("BEGIN QCSubDomain::printHeadOfInternalStructure " << this->overlapMap[NB_OVRLP_IDX] );
00837 std::cout <<" ############################################################################################" <<std::endl;
00838 std::cout <<" Head Of internal structure of subdomain "<< domainId << " numGlob : " << globalDomainId <<std::endl
00839 << " nbAtomCoreShell1 "<< this->nbAtomCoreShell1 <<std::endl
00840 << " nbAOCoreShell1 "<< this->nbAOCoreShell1 <<std::endl
00841 << " nbAtomShell2 "<< this->nbAtomShell2 <<std::endl
00842 << " nbAOShell2 "<< this->nbAOShell2 <<std::endl
00843 <<std::endl ;
00844 for ( n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00845 offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00846 std::cout << "SD num "<< n << " offset : " <<offset << " HEAD : " << std::endl;
00847 std::cout<< " SBD_ID_IDX : " << this->overlapMap[offset + SBD_ID_IDX ] << " " << std::endl;
00848 std::cout<< " NB_SHARED_ATOM : " << this->overlapMap[offset + NB_SHARED_ATOM_IDX ] << " " << std::endl;
00849 std::cout<< " NB_ATOM_SHELL1 : " << this->overlapMap[offset + NB_ATOM_SHELL1_IDX ] << " " << std::endl;
00850 std::cout<< " NB_ATOM_SHELL2 : " << this->overlapMap[offset + NB_ATOM_SHELL2_IDX ] << " " << std::endl;
00851 std::cout<< " NB_SHARED_AO : " << this->overlapMap[offset + NB_SHARED_AO_IDX ] << " " << std::endl;
00852 std::cout<< " NB_AO_SHELL1 : " << this->overlapMap[offset + NB_AO_SHELL1_IDX ] << " " << std::endl;
00853 std::cout<< " NB_AO_SHELL2 : " << this->overlapMap[offset + NB_AO_SHELL2_IDX ] << " " << std::endl;
00854 std::cout << std::endl;
00855 }
00856 std::cout <<" ############################################################################################" <<std::endl;
00857 QC_TRACE("END QCSubDomain::printHeadOfInternalStructure ");
00858 }
00862 void QCSubDomain::writeSubDomainHeader(std::ofstream & out){
00863 QC_TRACE("BEGIN QCSubDomain::writeSubDomainHeader ");
00864 int offset =0 ;
00865 int nbAtom , nbAO ;
00866
00867 nbAtom = this->nbAtomCoreShell1 + this->nbAtomShell2 ;
00868 nbAO = this->nbAOCoreShell1 + this->nbAOShell2 ;
00869 out << " "<< this->domainId << " " << nbAtom <<" " << nbAO << " "
00870 << this->overlapMap[NB_OVRLP_IDX] << " " ;
00871
00872
00873 for ( int n = 0 ; n < this->overlapMap[NB_OVRLP_IDX] ; ++n){
00874 offset = this->overlapMap[OVRLP_MAP_HEAD + n] ;
00875 out << this->overlapMap[offset + SBD_ID_IDX] <<" "<< this->overlapMap[offset + NB_SHARED_ATOM_IDX] << " ";
00876 }
00877 out << std::endl;
00878 QC_TRACE("END QCSubDomain::writeSubDomainHeader ");
00879 }
00880 void QCSubDomain::fillAtomInStructure(std::vector<std::vector<QCAtomIn> >& atomInZone){
00881
00882 QC_TRACE("BEGIN QCSubDomain::fillAtomInStructure ");
00883 int n ;
00884 QCAtomIn current;
00885
00886 current.numDomain = this->domainId ;
00887 for ( n = 0 ; n < this->getNbAtoms() ; ++n){
00888 current.localNum = n;
00889 current.typeZone = this->zone[n];
00890
00891 (atomInZone[this->globalIdx[n]]).push_back(current);
00892 }
00893 QC_TRACE("END QCSubDomain::fillAtomInStructure ");
00894 }
00895
00899 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::setSystemNumbers);
00900
00901 QCPARAMETER_METH_EXPL_INST_PARAM(void QCSubDomain::fillDensityArray,
00902 TWO_PARAMS(QCFloat * QCRestrict &, bool));
00903
00904 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::initWeightedDensity);
00905
00906 QCPARAMETER_METH_EXPL_INST(void QCSubDomain::computeBFactor);
00907
00908 QCPARAMETER_METH_EXPL_INST_PARAM(QCFloat QCSubDomain::traceProductCoreShell1,
00909 TWO_PARAMS(QCSymMatrix& QCRestrict, QCSymMatrix& QCRestrict));