00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <iomanip>
00021
00022 #include "config.h"
00023 #include "QCCommon.hpp"
00024 #include "QCFockGtr.hpp"
00025 #include "QCSymMatrix.hpp"
00026 #include "QCMndo.hpp"
00027 #include "QCAm1.hpp"
00028 #include "QCPm3.hpp"
00029 #include "QCMndoParam.hpp"
00030 #include "QCAm1Param.hpp"
00031 #include "QCPm3Param.hpp"
00032 #include "QCGlobalSystem.hpp"
00033 #include "QCMDSystem.hpp"
00034 #include "QCDistMDSystem.hpp"
00035 #include "QCSCFAlgo.hpp"
00036 #include "QCDCAlgo.hpp"
00037 #include "QCManager.hpp"
00038 #include "QCIntgReader.hpp"
00039
00040
00041
00045 template <class TPMatrix>
00046 QCFockGtr<TPMatrix>::QCFockGtr (void) :
00047 QCMatElemGtr<TPMatrix>()
00048 {}
00049
00050
00051
00055 template <class TPMatrix>
00056 QCFockGtr<TPMatrix>::QCFockGtr (int dimMatrix) :
00057 QCMatElemGtr<TPMatrix>(dimMatrix)
00058 {}
00059
00060
00061
00065 template <class TPMatrix>
00066 QCFockGtr<TPMatrix>::~QCFockGtr (void) {
00067 }
00068
00069
00070
00071
00079 template <class TPMatrix>
00080 template <class TPManager>
00081 void QCFockGtr<TPMatrix>::completeElems (TPManager& manager,
00082 typename TPManager::TSystem::QCIterator& workingSystem,
00083 const bool& isItfirstCall) {
00084
00085 const QCGeneralData& data = manager.getGeneralData();
00086
00087 QCIntgAcquisitionMethod intgAcquisitionMethod = data.getIntgAcquisitionMethod();
00088
00089
00090
00091
00092 QCSystem& QCRestrict currentSystem = workingSystem->getSystem();
00093 QCModelMatrices& QCRestrict currentMatrices = workingSystem->getMatrices();
00094
00095 const typename TPManager::TParam *QCRestrict parameters = manager.getParameters();
00096
00097 typename TPManager::TModel& QCRestrict model = manager.getModel();
00098
00099 QCPoint3D coordsA;
00100 QCPoint3D coordsB;
00101
00102 const typename TPManager::TParam *QCRestrict parameterA;
00103 const typename TPManager::TParam *QCRestrict parameterB;
00104
00105 QCFloat rAB;
00106
00107 QCSymMatrix& QCRestrict densityMatrix = currentMatrices.getDensityP().getMatrix();
00108
00109
00110 QCRepInterElec& QCRestrict repInterElecIntegrals = currentMatrices.getRepInterElecIntegrals();
00111
00112 QCMatrix& QCRestrict repIntegralsMatrix = repInterElecIntegrals.getMatrix();
00113
00114
00115
00116 QCSymMatrix& QCRestrict fockFAA = model.getSpWorkingAA();
00117 QCMatrix& QCRestrict fockFAB = model.getSpWorkingAB();
00118 QCSymMatrix& QCRestrict fockFBB = model.getSpWorkingBB();
00119
00120
00121 QCSymMatrix& QCRestrict densityPAA = model.getSpWorkingAA2();
00122 QCSymMatrix& QCRestrict densityPBB = model.getSpWorkingBB2();
00123 QCMatrix& QCRestrict densityPAB = model.getSpWorkingAB2();
00124
00125
00126
00127
00128 int firstAOofA, firstAOofB;
00129
00130 register QCFloat interVal, interVal1, interVal2;
00131 int i, j, mu, nu, lda, sma, muNu, ldaSma;
00132
00133
00134
00135 for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00136 currentSystem.getPointAt(i, coordsA);
00137 parameterA = ¶meters[currentSystem.getParamIndexAt(i)];
00138 firstAOofA = currentSystem.getFirstAOAt(i);
00139
00140
00141
00142 this->matrix.extractSubTriangle(fockFAA, firstAOofA, parameterA->getNbAO());
00143
00144
00145 densityMatrix.extractSubTriangle(densityPAA, firstAOofA, parameterA->getNbAO());
00146
00147
00148
00149 for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00150
00151 interVal = QC_ZERO;
00152 for (nu = 0; nu < parameterA->getNbAO(); ++nu) {
00153 interVal += densityPAA[nu][nu] * (parameterA->getG(mu, nu) - QC_HALF * parameterA->getH(mu, nu) );
00154 }
00155 fockFAA[mu][mu] += interVal;
00156 for (nu = 0; nu < mu; ++nu) {
00157 fockFAA[mu][nu] += QC_HALF * densityPAA[mu][nu] *
00158 ( QC_THREE * parameterA->getH(mu, nu) - parameterA->getG(mu, nu) );
00159 }
00160 }
00161
00162
00163
00164 for (j = 0; j < i; ++j) {
00165
00166 currentSystem.getPointAt(j, coordsB);
00167 parameterB = ¶meters[currentSystem.getParamIndexAt(j)];
00168 firstAOofB = currentSystem.getFirstAOAt(j);
00169 rAB = currentSystem.distance(i,j);
00170
00171
00172 this->matrix.extractSubMatrix(fockFAB,firstAOofA, parameterA->getNbAO(),
00173 firstAOofB, parameterB->getNbAO() );
00174 this->matrix.extractSubTriangle(fockFBB,firstAOofB, parameterB->getNbAO() );
00175
00176
00177 densityMatrix.extractSubMatrix(densityPAB, firstAOofA, parameterA->getNbAO(),
00178 firstAOofB, parameterB->getNbAO() );
00179 densityMatrix.extractSubTriangle(densityPBB, firstAOofB, parameterB->getNbAO() );
00180
00181
00182
00183 repInterElecIntegrals.computeElems(
00184
00185 parameterA,parameterB,coordsA,coordsB, rAB,
00186
00187
00188
00189 false, false, false, intgAcquisitionMethod, isItfirstCall);
00190
00191
00192
00193 for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00194 for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00195 interVal = QC_ZERO;
00196 muNu = mu * parameterA->getNbAO();
00197 for (nu = 0; nu < parameterA->getNbAO(); ++nu, ++muNu) {
00198 ldaSma = lda * parameterB->getNbAO();
00199 for (sma = 0; sma < parameterB->getNbAO(); ++sma, ++ldaSma) {
00200
00201
00202 interVal += densityPAB[nu][sma] * repIntegralsMatrix[muNu][ldaSma];
00203 }
00204 }
00205 fockFAB[mu][lda] -= QC_HALF * interVal;
00206 }
00207 }
00208
00209
00210
00211
00212
00213 if ( (data.getPartitionType() == QC_STANDARD_PART && workingSystem->getZone(j) != QC_SHELL2) ||
00214 (data.getPartitionType() == QC_DIXON_PART && workingSystem->getZone(j) == QC_CORE) ) {
00215 for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00216 muNu = mu * parameterA->getNbAO();
00217 for (nu = 0; nu <= mu; ++nu, ++muNu) {
00218 interVal1 = QC_ZERO;
00219 interVal2 = QC_ZERO;
00220 for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00221 ldaSma = lda * parameterB->getNbAO();
00222 for (sma = 0; sma < lda; ++sma, ++ldaSma) {
00223
00224
00225 interVal1 += densityPBB[lda][sma] * repIntegralsMatrix[muNu][ldaSma];
00226 }
00227
00228
00229 interVal2 += densityPBB[lda][sma] * repIntegralsMatrix[muNu][ldaSma];
00230 }
00231
00232
00233 fockFAA[mu][nu] += QC_TWO * interVal1;
00234 fockFAA[mu][nu] += interVal2;
00235 }
00236 }
00237 }
00238
00239
00240
00241
00242 if ( (data.getPartitionType() == QC_STANDARD_PART && workingSystem->getZone(i) != QC_SHELL2) ||
00243 (data.getPartitionType() == QC_DIXON_PART && workingSystem->getZone(i) == QC_CORE) ) {
00244
00245 for (lda = 0; lda < parameterB->getNbAO(); ++lda) {
00246 ldaSma = lda * parameterB->getNbAO();
00247 for (sma = 0; sma <= lda; ++sma, ++ldaSma) {
00248 interVal1 = QC_ZERO;
00249 interVal2 = QC_ZERO;
00250 for (mu = 0; mu < parameterA->getNbAO(); ++mu) {
00251 muNu = mu * parameterA->getNbAO();
00252 for (nu = 0; nu < mu; ++nu, ++muNu) {
00253
00254
00255 interVal1 += densityPAA[mu][nu] * repIntegralsMatrix[muNu][ldaSma];
00256 }
00257
00258
00259
00260 interVal2 += densityPAA[mu][nu] * repIntegralsMatrix[muNu][ldaSma];
00261 }
00262
00263
00264 fockFBB[lda][sma] += QC_TWO * interVal1;
00265 fockFBB[lda][sma] += interVal2;
00266 }
00267 }
00268 }
00269
00270 this->matrix.insertSubMatrix(fockFAB, firstAOofA, parameterA->getNbAO(),
00271 firstAOofB, parameterB->getNbAO() );
00272
00273
00274 this->matrix.insertSubTriangle(fockFBB, firstAOofB, parameterB->getNbAO() );
00275 }
00276
00277 this->matrix.insertSubTriangle(fockFAA, firstAOofA, parameterA->getNbAO() );
00278 }
00279
00280 if (!manager.isDCComputation() && intgAcquisitionMethod == QC_INDIRECT_STORAGE && !isItfirstCall) {
00281 repInterElecIntegrals.getIntgReader()->rewindFile();
00282 }
00283 #ifdef VERBOSE_FOCK_BEFORE
00284 if (isItfirstCall) {
00285 ostringstream osstr;
00286 int id = workingSystem->getId();
00287 osstr << "fock_mine_" << id << "_before";
00288 this->matrix.printInFile(osstr.str().c_str());
00289 }
00290 #endif
00291 }
00295 template <class TPMatrix>
00296 template <class TPParam>
00297 void QCFockGtr<TPMatrix>::completeElems (const TPParam * QCRestrict params,
00298 const QCSubDomain& QCRestrict workingDomain,
00299 const QCSubDomain& QCRestrict remoteDomain,
00300 const QCFloat * QCRestrict remoteDensity,
00301 QCRepInterElec& QCRestrict repInterElecIntegrals,
00302 QCSymMatrix& QCRestrict fockFAA,
00303 const QCGeneralData& data,
00304 int step,
00305 int nbDomains,
00306 QCIntgAcquisitionMethod intgAcquisitionMethod,
00307 bool isFirstCall)
00308 {
00309
00310
00311 QCMatrix& QCRestrict repIntegralsMatrix = repInterElecIntegrals.getMatrix();
00312
00313
00314 const TPParam * QCRestrict paramA;
00315 const TPParam * QCRestrict paramB;
00316
00317
00318 QCFloat rAB;
00319
00320
00321 QCPoint3D coordsA;
00322 QCPoint3D coordsB;
00323
00324
00325 int firstAOofA;
00326
00327
00328 register QCFloat interVal1, interVal2;
00329
00330
00331 int i, j, k, mu, nu, lda, sma, muNu, ldaSma;
00332
00333
00334 int blockIdx, firstBlockIdx;
00335
00336
00337 for (i = 0; i < workingDomain.getNbAtoms(); ++i) {
00338
00339 paramA = ¶ms[workingDomain.getParamIndexAt(i)];
00340 workingDomain.getPointAt(i, coordsA);
00341 firstAOofA = workingDomain.getFirstAOAt(i);
00342
00343 this->matrix.extractSubTriangle(fockFAA, firstAOofA, paramA->getNbAO());
00344
00345
00346 for (j = 0; j < remoteDomain.getNbAtoms(); ++j) {
00347
00348
00349 bool contributed = false;
00350
00351 if (data.getPartitionType() == QC_STANDARD_PART) {
00352
00353 if (remoteDomain.getZone(j) == QC_SHELL2) {
00354 contributed = true;
00355
00356 } else if (remoteDomain.getZone(j) == QC_SHELL1) {
00357
00358 for (k=0; k <= step; ++k) {
00359 int contId = (remoteDomain.getId() + k + 1) % nbDomains;
00360 int remoteIdx = remoteDomain.getRemoteIdx(contId, j);
00361
00362 if (remoteIdx >= 0) {
00363 contributed = true;
00364 break;
00365 }
00366 }
00367 }
00368
00369 }
00370 else if ( data.getPartitionType() == QC_DIXON_PART ) {
00371 if (remoteDomain.getZone(j) != QC_CORE) {
00372 contributed = true;
00373 }
00374 }
00375 if (!contributed &&
00376 (workingDomain.getRemoteIdx(remoteDomain.getId(), i) != j)) {
00377
00378 paramB = ¶ms[remoteDomain.getParamIndexAt(j)];
00379 remoteDomain.getPointAt(j, coordsB);
00380
00381
00382 rAB = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00383 QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00384 QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00385
00386
00387 repInterElecIntegrals.computeElems(paramA,paramB, coordsA,coordsB, rAB,
00388 false, false,false, intgAcquisitionMethod,
00389 isFirstCall);
00390
00391 firstBlockIdx = remoteDomain.getDensityBlockIdx(j);
00392
00393 for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00394 muNu = mu * paramA->getNbAO();
00395 for (nu = 0; nu <= mu; ++nu, ++muNu) {
00396 blockIdx = firstBlockIdx;
00397 interVal1 = QC_ZERO;
00398 interVal2 = QC_ZERO;
00399
00400 for (lda = 0; lda < paramB->getNbAO(); ++lda) {
00401 ldaSma = lda * paramB->getNbAO();
00402 for (sma = 0; sma < lda; ++sma, ++ldaSma, ++blockIdx) {
00403
00404
00405
00406 interVal1 += remoteDensity[blockIdx] * repIntegralsMatrix[muNu][ldaSma];
00407 }
00408
00409
00410
00411 interVal2 += remoteDensity[blockIdx] * repIntegralsMatrix[muNu][ldaSma];
00412 ++blockIdx;
00413 }
00414
00415
00416 fockFAA[mu][nu] += QC_TWO * interVal1;
00417 fockFAA[mu][nu] += interVal2;
00418 }
00419 }
00420 }
00421 }
00422
00423
00424 this->matrix.insertSubTriangle(fockFAA, firstAOofA, paramA->getNbAO());
00425 }
00426 }
00431 template <class TPMatrix>
00432 template <class TPIterator>
00433 void
00434 QCFockGtr<TPMatrix>::levelShifting (TPIterator& workingSystem, const QCFloat& doubleGapShift) {
00435
00436
00437 QCSymMatrix& QCRestrict densityMatrix = workingSystem->getMatrices().getDensityP().getMatrix();
00438
00439
00440 QCFloat gapShift = QC_HALF * doubleGapShift;
00441 int i, j;
00442
00443
00444
00445 for (i = 0; i < densityMatrix.getDim() ; ++i) {
00446 for (j = 0; j <= i; ++j) {
00447
00448 this->matrix[i][j] += gapShift * densityMatrix[i][j];
00449 }
00450
00451 this->matrix[i][i] -= doubleGapShift;
00452
00453
00454
00455
00456 }
00457 }
00458
00465 template <class TPMatrix>
00466 template <class TPManager>
00467 void
00468 QCFockGtr<TPMatrix>::diagonalize (TPManager& manager,
00469 typename TPManager::TSystem::QCIterator& workingSystem,
00470 const int& nbIter) {
00471
00472
00473 QCMemory& memory = manager.getMemory();
00474
00475
00476 QCModelMatrices& QCRestrict matrices =
00477 workingSystem->getMatrices();
00478
00479
00480 QCMatrix& transEigenVectMatrix = matrices.getEigenVectCt().getMatrix();
00481
00482
00483 QCFloat * QCRestrict eigenVal = matrices.getEigenVal();
00484
00485
00486 int * QCRestrict degener = matrices.getDegener();
00487 const QCFloat toleranceFactor = 1e-12;
00488
00489
00490 QCDiagoAlgorithm diagoAlgorithm =
00491 manager.getGeneralData().getDiagoAlgorithm();
00492
00493 #ifdef USE_LAPACK
00494
00495
00496 if (diagoAlgorithm == QC_DC_DIAGO) {
00497
00498
00499
00500
00501
00502
00503
00504 this->matrix.copyIn(transEigenVectMatrix);
00505 }
00506 #endif // USE_LAPACK
00507
00508
00509 if (!this->matrix.diagonalize(memory, diagoAlgorithm, toleranceFactor, eigenVal,
00510 degener,transEigenVectMatrix.getElems(), nbIter)) {
00511 cerr << __FILE__ << "diagonalise: diagonalisation aborted"
00512 << endl;
00513 exit(EXIT_FAILURE);
00514 }
00515 }
00516
00517
00518
00519
00520
00521
00522
00526 template class QCFockGtr<QCSymMatrix>;
00527
00528 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::completeElems,
00529 const bool&);
00530
00531 QCPARAMETER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::completeElems,
00532 TEN_PARAMS(const QCSubDomain& QCRestrict,
00533 const QCSubDomain& QCRestrict,
00534 const QCFloat * QCRestrict,
00535 QCRepInterElec& QCRestrict,
00536 QCSymMatrix& QCRestrict,
00537 const QCGeneralData&,
00538 int,
00539 int,
00540 QCIntgAcquisitionMethod,
00541 bool));
00542
00543 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::diagonalize,
00544 const int&);
00545
00546 QCITERATOR_METH_EXPL_INST_PARAM(void QCFockGtr<QCSymMatrix>::levelShifting,
00547 const QCFloat&);
00548