00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "QCCommon.hpp"
00022 #include "QCTable.hpp"
00023 #include "QCDensityGtr.hpp"
00024 #include "QCSymMatrix.hpp"
00025 #include "QCSystem.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 "QCSubDomain.hpp"
00039
00040 #if defined (QC_TRACE_IN_FILE)
00041 #include "QCTrace.hpp"
00042 #endif
00043
00047 template <class TPMatrix>
00048 QCDensityGtr<TPMatrix>::QCDensityGtr (void) :
00049 QCMatElemGtr<TPMatrix>()
00050 {}
00051
00052
00053
00054
00058 template <class TPMatrix>
00059 QCDensityGtr<TPMatrix>::QCDensityGtr (int dimMatrix) :
00060 QCMatElemGtr<TPMatrix>(dimMatrix)
00061 {}
00062
00063
00064
00065
00069 template <class TPMatrix>
00070 QCDensityGtr<TPMatrix>::~QCDensityGtr (void) {
00071 }
00072
00073 template <class TPMatrix>
00074 template <class TPManager>
00075 void
00076 QCDensityGtr<TPMatrix>::initDiag (TPManager& manager,
00077 typename TPManager::TSystem::QCIterator& iter) {
00078 QC_TRACE_INIT("BEGIN QCDensityGtr<TPMatrix>::initDiag") ;
00079
00080
00081 QCSystem& QCRestrict currentSystem = iter->getSystem();
00082 const QCSystem& QCRestrict rootSystem = iter->getRootSystem();
00083
00084
00085 const typename TPManager::TParam *QCRestrict parameters =
00086 manager.getParameters();
00087
00088
00089 const typename TPManager::TParam *QCRestrict parameterA;
00090
00091
00092
00093 int firstAOofA, firstAOofSucA;
00094
00095 QCFloat diagElem, chargeForDiagElem;
00096
00097 int i, mu;
00098
00099 chargeForDiagElem = ( static_cast<QCFloat>(rootSystem.getSystemCharge()) /
00100 rootSystem.getNbAtomicOrbitals() );
00101 this->matrix.initClean();
00102
00103 for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00104
00105 parameterA = ¶meters[currentSystem.getParamIndexAt(i)];
00106 firstAOofA = currentSystem.getFirstAOAt(i);
00107 firstAOofSucA = firstAOofA + parameterA->getNbAO();
00108
00109
00110
00111 diagElem = ( ( static_cast<QCFloat>(parameterA->getNbValenceElec()) /
00112 parameterA->getNbAO() )
00113
00114 - chargeForDiagElem );
00115 for (mu = firstAOofA; mu < firstAOofSucA; ++mu) {
00116 this->matrix[mu][mu] = diagElem;
00117 }
00118 }
00119
00120 QC_TRACE_INIT("END QCDensityGtr<TPMatrix>::initDiag") ;
00121 }
00125 template <class TPMatrix>
00126 template <class TPMatrix2>
00127 void
00128 QCDensityGtr<TPMatrix>::computeElems (QCMemory& QCRestrict memory,
00129 const int nbElecOfSystem,
00130 QCMatElemGtr<TPMatrix2>& QCRestrict eigenVectC) {
00131
00132
00133 TPMatrix2& QCRestrict transEigenVectMatrix = eigenVectC.getMatrix();
00134
00135 const int nbFilledOM = nbElecOfSystem / 2;
00136 const int dim = this->matrix.getDim();
00137
00138 int i;
00139 QCFloat* QCRestrict workingMatrix =
00140 memory.template takeTmpMem<QCFloat>(dim*dim);
00141
00148
00149
00150 QCCBlas(dsyrk)(
00151 #ifdef USE_OPTIMIZED_CBLAS
00152 CblasRowMajor,
00153 CblasLower, CblasTrans,
00154 #else // USE_OPTIMIZED_CBLAS
00155 CblasUpper, CblasNoTrans,
00156 #endif // USE_OPTIMIZED_CBLAS
00157 dim, nbFilledOM,
00158 QC_TWO, transEigenVectMatrix.getElems(), dim,
00159 QC_ZERO, workingMatrix, dim);
00160
00161
00162
00163 for (i = 0; i < dim; ++i) {
00164 memcpy(this->matrix[i], &(workingMatrix[i*dim]), (i+1) * sizeof(QCFloat) );
00165 }
00166
00167 memory.template giveBackTmpMem<QCFloat>(workingMatrix);
00168 }
00172 template <class TPMatrix>
00173 template <class TPParam>
00174 void
00175 QCDensityGtr<TPMatrix>::mergeDensity (const TPParam * QCRestrict parameters,
00176 QCSubDomain& QCRestrict workingDomain,
00177 QCSubDomain& QCRestrict remoteDomain,
00178 QCSymMatrix& QCRestrict remWeightedDensity,
00179 int remoteMapIdx,
00180 QCMatrix& QCRestrict wDensityAB,
00181 QCMatrix& QCRestrict interPAB,
00182 QCSymMatrix& QCRestrict wDensityAA) {
00183
00184 QC_TRACE("BEGIN QCDensityGtr<TPMatrix>::mergeDensity ");
00185
00186
00187
00188 const TPParam * QCRestrict paramA;
00189 const TPParam * QCRestrict paramB;
00190
00191
00192 int firstAOofA, firstAOofB;
00193 int nbAOofA, nbAOofB;
00194 int localFirstAOHigh, localFirstAOLow;
00195 int localNbAOLow, localNbAOHigh;
00196 int remoteNbAOLow, remoteNbAOHigh;
00197 int remoteLow, remoteHigh;
00198 bool atomAisHighRemote, atomAisHighLocal;
00199
00200
00201 int localMu, localNu, localLda;
00202 int mu, nu, lda, a, b;
00203
00204 for (a = 0; a < workingDomain.getNbSharedAtoms(remoteMapIdx); ++a) {
00205
00206 int localA = workingDomain.getIndLocalIdx(remoteMapIdx, a);
00207 int remoteA = workingDomain.getIndRemoteIdx(remoteMapIdx, a);
00208 paramA = ¶meters[workingDomain.getParamIndexAt(localA)];
00209 firstAOofA = workingDomain.getFirstAOAt(localA);
00210 nbAOofA = paramA->getNbAO();
00211
00212 #ifdef DEBUG_MERGE_DENSITY
00213 cout << endl << "theDiDkCellPtrIA->localNumInThis = " << localA << " pour " << table[paramA->getIdentifier()] << endl;
00214 cout << "theDiDkCellPtrIA->localNumInOther = " << remoteA << endl;
00215 #endif
00216
00217 for (b = 0; b < a; ++b) {
00218
00219 int localB = workingDomain.getIndLocalIdx(remoteMapIdx, b);
00220 int remoteB = workingDomain.getIndRemoteIdx(remoteMapIdx, b);
00221 paramB = ¶meters[workingDomain.getParamIndexAt(localB)];
00222 firstAOofB = workingDomain.getFirstAOAt(localB);
00223 nbAOofB = paramB->getNbAO();
00224
00225 if (firstAOofA > firstAOofB) {
00226
00227 localFirstAOHigh = firstAOofA;
00228 localFirstAOLow = firstAOofB;
00229 localNbAOHigh = nbAOofA;
00230 localNbAOLow = nbAOofB;
00231 atomAisHighLocal = true;
00232
00233 } else {
00234
00235 localFirstAOHigh = firstAOofB;
00236 localFirstAOLow = firstAOofA;
00237 localNbAOHigh = nbAOofB;
00238 localNbAOLow = nbAOofA;
00239 atomAisHighLocal = false;
00240
00241 }
00242
00243 if (remoteA > remoteB) {
00244
00245 remoteHigh = remoteA;
00246 remoteLow = remoteB;
00247 remoteNbAOHigh = nbAOofA;
00248 remoteNbAOLow = nbAOofB;
00249 atomAisHighRemote = true;
00250
00251 } else {
00252
00253 remoteHigh = remoteB;
00254 remoteLow = remoteA;
00255 remoteNbAOHigh = nbAOofB;
00256 remoteNbAOLow = nbAOofA;
00257 atomAisHighRemote = false;
00258
00259 }
00260
00261
00262 remWeightedDensity.extractSubMatrix(wDensityAB,
00263 remoteDomain.getFirstAOAt(remoteHigh),
00264 remoteNbAOHigh,
00265 remoteDomain.getFirstAOAt(remoteLow),
00266 remoteNbAOLow);
00267
00268
00269
00270 if ( (atomAisHighLocal && !atomAisHighRemote) ||
00271 (!atomAisHighLocal && atomAisHighRemote) ) {
00272
00273 interPAB.transpose(wDensityAB,
00274 remoteNbAOHigh,
00275 remoteNbAOLow);
00276 interPAB.swap(wDensityAB);
00277 }
00278
00279
00280 for (mu = 0, localMu = localFirstAOHigh; mu < localNbAOHigh; ++mu, ++localMu) {
00281 for (lda = 0, localLda = localFirstAOLow; lda < localNbAOLow; ++lda, ++localLda) {
00282
00283 #ifdef DEBUG_MERGE_DENSITY
00284 cout << "avant : mergedDensityMatrix[ " << localMu << "][" << localLda << "] = "
00285 << this->matrix[localMu][localLda]
00286 << " et weightedDensityPABInOther[" << mu << "][ " << lda << "] = "
00287 << wDensityAB[mu][lda] << endl;
00288 #endif
00289
00290 this->matrix[localMu][localLda] += wDensityAB[mu][lda];
00291
00292 #ifdef DEBUG_MERGE_DENSITY
00293 cout << "apres : mergedDensityMatrix[ " << localMu << "][" << localLda << "] = "
00294 << this->matrix[localMu][localLda]
00295 << " et weightedDensityPABInOther[" << mu << "][ " << lda << "] = "
00296 << wDensityAB[mu][lda] << endl;
00297 #endif
00298
00299 }
00300 }
00301 }
00302
00303
00304 remWeightedDensity.extractSubTriangle(wDensityAA, remoteDomain.getFirstAOAt(remoteA), nbAOofA);
00305
00306 for (mu = 0, localMu = firstAOofA; mu < nbAOofA; ++mu, ++localMu) {
00307 for (nu = 0, localNu = firstAOofA; nu <= mu; ++nu, ++localNu) {
00308 this->matrix[localMu][localNu] += wDensityAA[mu][nu];
00309 }
00310 }
00311 }
00312 QC_TRACE("END QCDensityGtr<TPMatrix>::mergeDensity ");
00313 }
00317 template <class TPMatrix>
00318 template <class TPMatrix2>
00319 void
00320 QCDensityGtr<TPMatrix>::computeElemsDC (QCMemory& QCRestrict memory,
00321 const int nbFilledOM,
00322 const int nbOccupiedOM,
00323 const int nbAOCoreShell1,
00324 const QCFloat * QCRestrict orbitalOccupN,
00325 QCMatElemGtr<TPMatrix2>& QCRestrict eigenVectCt) {
00326
00327
00328
00329 TPMatrix2& QCRestrict transEigenVectMatrix = eigenVectCt.getMatrix();
00330
00331 const int dim = this->matrix.getDim();
00332 int i;
00333 QCFloat* QCRestrict workingMatrix = memory.template takeTmpMem<QCFloat>(dim*dim);
00334
00335 #define QC_CUT_DENSITY
00336
00343 for (i = nbFilledOM; i < nbOccupiedOM; ++i) {
00344
00345 #ifdef QC_CUT_DENSITY
00346 QCCBlas(dscal)(nbAOCoreShell1, sqrt(orbitalOccupN[i] / QC_TWO),
00347 transEigenVectMatrix[i], 1);
00348 #else // QC_CUT_DENSITY
00349 QCCBlas(dscal)(dim, sqrt(orbitalOccupN[i] / QC_TWO),
00350 transEigenVectMatrix[i], 1);
00351 #endif // QC_CUT_DENSITY
00352
00353 }
00360
00361
00362 QCCBlas(dsyrk)(
00363 #ifdef USE_OPTIMIZED_CBLAS
00364 CblasRowMajor,
00365 CblasLower, CblasTrans,
00366 #else // USE_OPTIMIZED_CBLAS
00367 CblasUpper, CblasNoTrans,
00368 #endif // USE_OPTIMIZED_CBLAS
00369
00370
00371 #ifdef QC_CUT_DENSITY
00372 nbAOCoreShell1, nbOccupiedOM,
00373 #else // QC_CUT_DENSITY
00374 dim, nbOccupiedOM,
00375 #endif // QC_CUT_DENSITY
00376 QC_TWO, transEigenVectMatrix.getElems(), dim,
00377 QC_ZERO, workingMatrix, dim);
00378
00379
00380
00381 #ifdef QC_CUT_DENSITY
00382 for (i = 0; i < nbAOCoreShell1; ++i) {
00383 #else // QC_CUT_DENSITY
00384 for (i = 0; i < dim; ++i) {
00385 #endif // QC_CUT_DENSITY
00386 memcpy(this->matrix[i], &(workingMatrix[i*dim]), (i+1) * sizeof(QCFloat) );
00387 }
00388 memory.template giveBackTmpMem<QCFloat>(workingMatrix);
00389 }
00390
00394 template class QCDensityGtr<QCSymMatrix>;
00395
00396 QCMANAGER_ITER_METH_EXPL_INST(void QCDensityGtr<QCSymMatrix>::initDiag);
00397
00398 template void
00399 QCDensityGtr<QCSymMatrix>::computeElems<QCMatrix>(QCMemory& QCRestrict ,
00400 const int ,
00401 QCMatElemGtr<QCMatrix>& QCRestrict);
00402
00403 template void
00404 QCDensityGtr<QCSymMatrix>::computeElemsDC<QCMatrix> (QCMemory& QCRestrict memory,
00405 const int nbFilledOM,
00406 const int nbOccupiedOM,
00407 const int nbAOCoreShell1,
00408 const QCFloat * QCRestrict orbitalOccupN,
00409 QCMatElemGtr<QCMatrix>& QCRestrict eigenVectCt);
00410
00411
00412 QCPARAMETER_METH_EXPL_INST_PARAM(void QCDensityGtr<QCSymMatrix>::mergeDensity,
00413 SEVEN_PARAMS(QCSubDomain& QCRestrict,
00414 QCSubDomain& QCRestrict,
00415 QCSymMatrix& QCRestrict,
00416 int,
00417 QCMatrix& QCRestrict,
00418 QCMatrix& QCRestrict,
00419 QCSymMatrix& QCRestrict));