00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include <sstream>
00020
00021 #include "QCCommon.hpp"
00022 #include "QCModelMatrices.hpp"
00023 #include "QCGlobalSystem.hpp"
00024 #include "QCMDSystem.hpp"
00025 #include "QCDistMDSystem.hpp"
00026 #include "QCMndo.hpp"
00027 #include "QCMndoParam.hpp"
00028 #include "QCAm1.hpp"
00029 #include "QCAm1Param.hpp"
00030 #include "QCPm3.hpp"
00031 #include "QCPm3Param.hpp"
00032 #include "QCSCFAlgo.hpp"
00033 #include "QCDCAlgo.hpp"
00034 #include "QCManager.hpp"
00035 #include "QCIntgWriter.hpp"
00036 #include "QCIntgReader.hpp"
00037
00038 static const char* IdempotFromEnv = getenv("IDEMPOT");
00039
00040 QCModelMatrices::QCModelMatrices (void) :
00041 overlapS(NULL),
00042 hamiltonH(),
00043 densityP(),
00044 densityPnm1(),
00045 fockF(),
00046 fockFTild(NULL),
00047 densityPTild(NULL),
00048 eigenVectCt(),
00049 eigenVal(NULL),
00050 degener(NULL),
00051 weightedDensityP(NULL),
00052 squareDensityP(NULL),
00053 threeIdMDensityP(NULL),
00054 repInterElecIntegrals(NULL),
00055 dRepInterElecIntgOnDCoord(NULL),
00056 dOverlapSOnDr(NULL)
00057 {
00058
00064 overlapS = new QCOverlapGtr<QCSymMatrix>();
00065
00066 }
00067
00071 QCModelMatrices::~QCModelMatrices (void) {
00072
00073
00074 if (overlapS) {
00075 delete overlapS; overlapS = NULL;
00076 }
00077 if (eigenVal) {
00078 delete [] eigenVal; eigenVal = NULL;
00079 }
00080 if (degener) {
00081 delete [] degener; degener = NULL;
00082 }
00083 if (weightedDensityP) {
00084 delete weightedDensityP; weightedDensityP = NULL;
00085 }
00086 if (squareDensityP) {
00087 delete squareDensityP; squareDensityP = NULL;
00088 }
00089 if (threeIdMDensityP) {
00090 delete threeIdMDensityP; threeIdMDensityP = NULL;
00091 }
00092
00093
00094 if (fockFTild) {
00095 delete fockFTild; fockFTild = NULL;
00096 }
00097 if (densityPTild) {
00098 delete densityPTild; densityPTild = NULL;
00099 }
00100
00101
00102 if (repInterElecIntegrals) {
00103 delete repInterElecIntegrals; repInterElecIntegrals = NULL;
00104 }
00105 if (dRepInterElecIntgOnDCoord) {
00106 delete dRepInterElecIntgOnDCoord; dRepInterElecIntgOnDCoord = NULL;
00107 }
00108 if (dOverlapSOnDr) {
00109 delete dOverlapSOnDr; dOverlapSOnDr = NULL;
00110 }
00111 }
00112
00116 template <class TPManager>
00117 void
00118 QCModelMatrices::allocateMatrices (TPManager& manager, int nbAOofSystem) {
00119
00120 typedef typename TPManager::TParam TParam;
00121
00122 const QCGeneralData& data = manager.getGeneralData();
00123
00124 bool useOptDump = ( data.getMainAlgorithm() == QC_OPTIMAL_DAMPING ||
00125 data.getMainAlgorithm() == QC_DC_OPTIMAL_DAMPING );
00126 bool isDCComputation = manager.isDCComputation();
00127 bool useDerivation = data.getComputationType() == QC_ENERGY_DERIV;
00128 int nbAOspOfAtom = ( TParam::NB_S_ORB + TParam::NB_P_ORB );
00129
00130
00131 overlapS->sizeMatrix(nbAOofSystem);
00132 hamiltonH.allocateMatrix(nbAOofSystem);
00133 densityP.allocateMatrix(nbAOofSystem);
00134
00135 #ifdef DEV_SCF_CLASS_NO_MEMORY
00136 densityPnm1.allocateMatrix(nbAOofSystem);
00137 #endif
00138 fockF.allocateMatrix(nbAOofSystem);
00139 eigenVectCt.allocateMatrix(nbAOofSystem);
00140
00141
00142 if (useOptDump) {
00143
00144 #ifdef DEV_SCF_CLASS_NO_MEMORY
00145
00146 fockFTild = new QCFockGtr<QCSymMatrix>(nbAOofSystem);
00147 densityPTild = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00148 #else
00149
00150
00151 fockFTild = new QCFockGtr<QCSymMatrix>();
00152 densityPTild = new QCDensityGtr<QCSymMatrix>();
00153 #endif
00154
00155 }
00156 eigenVal = new QCFloat [nbAOofSystem];
00157 degener = new int [nbAOofSystem];
00158
00159
00160 if (isDCComputation) {
00161
00162
00163
00164 weightedDensityP = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00165
00166 if (IdempotFromEnv) {
00167
00168 squareDensityP = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00169 threeIdMDensityP = new QCDensityGtr<QCSymMatrix>(nbAOofSystem);
00170 }
00171 }
00172
00173 repInterElecIntegrals = new QCRepInterElec (nbAOspOfAtom*nbAOspOfAtom);
00174
00175 if (useDerivation) {
00176
00177
00178 if (data.getDerivationType() == QC_NUMERICAL_DERIVATION) {
00179 dRepInterElecIntgOnDCoord = new QCRepInterElec(nbAOspOfAtom * nbAOspOfAtom);
00180 dOverlapSOnDr = new QCOverlapGtr<QCSymMatrix>();
00181 dOverlapSOnDr->sizeMatrix(nbAOofSystem);
00182 }
00183 }
00184 }
00185
00189 template <class TPManager>
00190 void
00191 QCModelMatrices::buildMatrices (TPManager& manager, typename TPManager::TSystem::QCIterator& workingSystem) {
00192
00193 QCIntgWriter *intgWriter = repInterElecIntegrals->getIntgWriter();
00194
00195
00196 const QCGeneralData& data = manager.getGeneralData();
00197 bool isDCComputation = manager.isDCComputation();
00198
00199
00200
00201 ostringstream idstrm;
00202 if (isDCComputation) {
00203 idstrm << workingSystem->getId();
00204 }
00205
00206 string integralsTmpFile = ( QCCommon::IntegralsTmpFile + "_" +
00207 idstrm.str() + "_" + OUTPUT_SUFFIX );
00208
00209
00210
00211
00212 if (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE) {
00217 if (!intgWriter) {
00218 intgWriter = repInterElecIntegrals->allocIntgWriter();
00219 repInterElecIntegrals->setIntgTmpFile(integralsTmpFile);
00220 }
00221 if (intgWriter) {
00222 intgWriter->openFile (integralsTmpFile.c_str(), "wb");
00223
00224 } else {
00225 cerr << __FILE__ << ": buildMatrices: error openning integrals file" << endl;
00226 exit(EXIT_FAILURE);
00227 }
00228 }
00229
00230 QCMemory& QCRestrict memory = manager.getMemory();
00231
00232
00233
00234 if (data.getInitDensityMatrix() == QC_DENSITY_DIAG) {
00235 densityP.initDiag(manager, workingSystem);}
00236
00237
00238
00239 overlapS->getMatrix().takeTmpMem(memory);
00240 overlapS->computeElems(manager, workingSystem);
00241
00242 memory.evaluateUsedMemory ();
00243
00244
00245
00246 hamiltonH.init();
00247 hamiltonH.computeElems(manager, workingSystem);
00248
00249 overlapS->getMatrix().giveBackTmpMem(memory);
00250
00251 if (!manager.isDCComputation()) {
00252
00253 fockF.copyMatrix(hamiltonH);
00254 fockF.completeElems(manager, workingSystem, true);
00255 }
00256 }
00257
00261 template <class TPManager>
00262 void QCModelMatrices::prepareIterations (TPManager& manager) {
00263
00264 QCIntgReader *intgReader = repInterElecIntegrals->getIntgReader();
00265 QCIntgWriter *intgWriter = repInterElecIntegrals->getIntgWriter();
00266
00267
00268 const QCGeneralData& data = manager.getGeneralData();
00269
00270
00271 QCMemory& memory = manager.getMemory();
00272
00273
00274 memory.updateBlockTmpMem();
00275
00276
00277
00278 if (data.getIntgAcquisitionMethod() == QC_INDIRECT_STORAGE) {
00279
00280
00281
00282 intgWriter->closeFile();
00283
00284
00285
00286
00287
00288
00289 if (!intgReader) {
00290 intgReader = repInterElecIntegrals->allocIntgReader();
00291 }
00292
00293 if (intgReader) {
00294 intgReader->openFile(repInterElecIntegrals->getIntgTmpFile().c_str(), "rb");
00295
00296 } else {
00297 cerr << __FILE__ << ": prepareIterations: error openning integrals file"
00298 << endl;
00299 exit(EXIT_FAILURE);
00300 }
00301 }
00302 }
00303
00307 QCMANAGER_METH_EXPL_INST_PARAM(void QCModelMatrices::allocateMatrices,
00308 int);
00309
00310 QCMANAGER_ITER_METH_EXPL_INST(void QCModelMatrices::buildMatrices);
00311
00312 QCMANAGER_METH_EXPL_INST(void QCModelMatrices::prepareIterations);