00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "QCCommon.hpp"
00023 #include "QCHamiltonGtr.hpp"
00024 #include "QCSymMatrix.hpp"
00025 #include "QCRepInterElec.hpp"
00026 #include "QCManager.hpp"
00027 #include "QCMndo.hpp"
00028 #include "QCAm1.hpp"
00029 #include "QCPm3.hpp"
00030 #include "QCMndoParam.hpp"
00031 #include "QCAm1Param.hpp"
00032 #include "QCPm3Param.hpp"
00033 #include "QCGlobalSystem.hpp"
00034 #include "QCMDSystem.hpp"
00035 #include "QCDistMDSystem.hpp"
00036 #include "QCSCFAlgo.hpp"
00037 #include "QCDCAlgo.hpp"
00038
00039
00040
00044 template <class TPMatrix>
00045 QCHamiltonGtr<TPMatrix>::QCHamiltonGtr (void) :
00046 QCMatElemGtr<TPMatrix>()
00047 {}
00048
00049
00050
00054 template <class TPMatrix>
00055 QCHamiltonGtr<TPMatrix>::~QCHamiltonGtr (void)
00056 {}
00057
00058
00059
00060
00064 template <class TPMatrix>
00065 template <class TPManager>
00066 void
00067 QCHamiltonGtr<TPMatrix>::computeElems (TPManager& manager,
00068 typename TPManager::TSystem::QCIterator& workingSystem) {
00069
00070
00071 typedef typename TPManager::TModel TModel;
00072 typedef typename TPManager::TModel::TParam TParam;
00073
00078 const QCSystem& QCRestrict currentSystem = workingSystem->getSystem();
00079 TModel& QCRestrict model = manager.getModel();
00080 const TParam * QCRestrict params = manager.getParameters();
00081
00082
00083 const QCGeneralData& data = manager.getGeneralData();
00084
00085
00086 const TParam * QCRestrict paramA;
00087 const TParam * QCRestrict paramB;
00088
00089
00090 QCFloat rAB;
00091
00092
00093 QCSymMatrix& QCRestrict overlapMatrix =
00094 workingSystem->getMatrices().getOverlapS()->getMatrix();
00095
00096
00097 QCRepInterElec& QCRestrict repInterElecIntegrals =
00098 workingSystem->getMatrices().getRepInterElecIntegrals();
00099
00100
00101 QCMatrix& QCRestrict repIntegralsMatrix =
00102 repInterElecIntegrals.getMatrix();
00103
00104
00105
00106 QCSymMatrix& QCRestrict hamiltonHAA = model.getSpWorkingAA();
00107 QCMatrix& QCRestrict hamiltonHAB = model.getSpWorkingAB();
00108
00109
00110 QCMatrix& QCRestrict overlapSAB = model.getSpdWorkingAB();
00111
00112
00113 const int SSRepIndex = 0;
00114
00115
00116
00117
00118 int firstAOofA, firstAOofB;
00119
00120
00121 QCPoint3D coordsA, coordsB;
00122
00123
00124 int i, j, mu, nu, lda, muNu;
00125
00126
00127 for (i = 0; i < currentSystem.getNbAtoms(); ++i) {
00128
00129 paramA = ¶ms[currentSystem.getParamIndexAt(i)];
00130 firstAOofA = currentSystem.getFirstAOAt(i);
00131 currentSystem.getPointAt(i, coordsA);
00132
00133
00134 hamiltonHAA.initClean();
00135
00136
00137
00138 for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00139
00140
00141 hamiltonHAA[mu][mu] = paramA->getU(mu);
00142 }
00143
00144
00145 for (j = 0; j < currentSystem.getNbAtoms(); ++j) {
00146 if (j != i) {
00147
00148 if ( (data.getPartitionType() == QC_STANDARD_PART &&
00149 workingSystem->getZone(j) != QC_SHELL2) ||
00150
00151 (data.getPartitionType() == QC_DIXON_PART &&
00152 workingSystem->getZone(j) == QC_CORE) ) {
00153
00154
00155
00156
00157
00158
00159
00160
00161 paramB = ¶ms[currentSystem.getParamIndexAt(j)];
00162 currentSystem.getPointAt(j, coordsB);
00163 rAB = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00164 QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00165 QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00166
00167 repInterElecIntegrals.computeElems(paramB,
00168 paramA,
00169 coordsB,
00170 coordsA,
00171 rAB,
00172 true,
00173 false);
00174
00175
00176
00177
00178 for (mu = 0, muNu = 0; mu < paramA->getNbAO(); ++mu) {
00179 muNu = mu * paramA->getNbAO();
00180 for (nu = 0; nu <= mu; ++nu, ++muNu) {
00181 hamiltonHAA[mu][nu] -=
00182 paramB->getCoreCharge() * repIntegralsMatrix[SSRepIndex][muNu];
00183 }
00184 }
00185 }
00186 }
00187 }
00188
00189
00190
00191
00192
00193 this->matrix.insertSubTriangle(hamiltonHAA,
00194 firstAOofA, paramA->getNbAO());
00195
00196
00197
00198
00199 for (j = 0; j < i; ++j) {
00200
00201 paramB = ¶ms[currentSystem.getParamIndexAt(j)];
00202 firstAOofB = currentSystem.getFirstAOAt(j);
00203
00204
00205 overlapMatrix.extractSubMatrix(overlapSAB,
00206 firstAOofA, paramA->getNbAO(),
00207 firstAOofB, paramB->getNbAO());
00208
00209
00210
00211 for (mu = 0; mu < paramA->getNbAO(); ++mu) {
00212 for (lda = 0; lda < paramB->getNbAO(); ++lda) {
00213
00214
00215 hamiltonHAB[mu][lda] =
00216 ( (paramA->getBeta(mu) + paramB->getBeta(lda) ) *
00217 QC_HALF * overlapSAB[mu][lda] );
00218 }
00219 }
00220
00221
00222 this->matrix.insertSubMatrix(hamiltonHAB,
00223 firstAOofA, paramA->getNbAO(),
00224 firstAOofB, paramB->getNbAO() );
00225
00226 }
00227 }
00228 #ifdef VERBOSE_HAMILTON_BEFORE
00229 ostringstream osstr;
00230 int id = workingSystem->getId();
00231 osstr << "hamilton_mine_";
00232 if (id >= 0) {
00233 osstr << id << "_before";
00234 }
00235 this->matrix.printInFile(osstr.str().c_str());
00236 #endif
00237 }
00238
00239
00240
00241
00242
00246 template <class TPMatrix>
00247 template <class TPParam>
00248 void
00249 QCHamiltonGtr<TPMatrix>::completeElems (const TPParam * QCRestrict params,
00250 const QCSubDomain& QCRestrict workingDomain,
00251 const QCSubDomain& QCRestrict remoteDomain,
00252 QCRepInterElec& QCRestrict repInterElecIntegrals,
00253 QCSymMatrix& QCRestrict hamiltonHAA,
00254 const QCGeneralData& data,
00255 int step,
00256 int nbDomains)
00257 {
00258
00259
00260
00261 QCMatrix& QCRestrict repIntegralsMatrix =
00262 repInterElecIntegrals.getMatrix();
00263
00264
00265 const TPParam * QCRestrict paramA;
00266 const TPParam * QCRestrict paramB;
00267
00268
00269 QCFloat rAB;
00270
00271
00272 const int SSRepIndex = 0;
00273
00274
00275 QCPoint3D coordsA;
00276 QCPoint3D coordsB;
00277
00278
00279 int i, j, k, mu, nu, muNu;
00280
00281
00282 for (i = 0; i < workingDomain.getNbAtoms(); ++i) {
00283
00284 paramA = ¶ms[workingDomain.getParamIndexAt(i)];
00285 workingDomain.getPointAt(i, coordsA);
00286
00287
00288
00289
00290
00291 hamiltonHAA.initClean();
00292
00293
00294 for (j = 0; j < remoteDomain.getNbAtoms(); ++j) {
00295
00296
00297 bool contributed = false;
00298
00299 if (data.getPartitionType() == QC_STANDARD_PART) {
00300
00301 if (remoteDomain.getZone(j) == QC_SHELL2) {
00302 contributed = true;
00303
00304 } else if (remoteDomain.getZone(j) == QC_SHELL1) {
00305
00306 for (k=0; k <= step; ++k) {
00307 int contId = (remoteDomain.getId() + k + 1) % nbDomains;
00308 int remoteIdx = remoteDomain.getRemoteIdx(contId, j);
00309
00310 if (remoteIdx >= 0) {
00311 contributed = true;
00312 break;
00313 }
00314 }
00315 }
00316
00317 } else if ( data.getPartitionType() == QC_DIXON_PART ) {
00318
00319 if (remoteDomain.getZone(j) != QC_CORE) {
00320 contributed = true;
00321 }
00322 }
00323
00324
00325 if (!contributed &&
00326 (workingDomain.getRemoteIdx(remoteDomain.getId(), i) != j)) {
00327
00328
00329
00330
00331
00332
00333
00334 paramB = ¶ms[remoteDomain.getParamIndexAt(j)];
00335 remoteDomain.getPointAt(j, coordsB);
00336 rAB = sqrt(QCPow<2>(coordsA[COORDX] - coordsB[COORDX]) +
00337 QCPow<2>(coordsA[COORDY] - coordsB[COORDY]) +
00338 QCPow<2>(coordsA[COORDZ] - coordsB[COORDZ]) );
00339
00340 repInterElecIntegrals.computeElems(paramB,
00341 paramA,
00342 coordsB,
00343 coordsA,
00344 rAB,
00345 true,
00346 false);
00347
00348
00349
00350 for (mu = 0, muNu = 0; mu < paramA->getNbAO(); ++mu) {
00351 muNu = mu * paramA->getNbAO();
00352 for (nu = 0; nu <= mu; ++nu, ++muNu) {
00353 hamiltonHAA[mu][nu] -=
00354 paramB->getCoreCharge() * repIntegralsMatrix[SSRepIndex][muNu];
00355 }
00356 }
00357 }
00358 }
00359
00360
00361
00362 this->matrix.completeSubTriangle(hamiltonHAA,
00363 workingDomain.getFirstAOAt(i),
00364 paramA->getNbAO());
00365 }
00366
00367 }
00368
00369
00370
00371
00372
00373
00374
00375
00379 template class QCHamiltonGtr<QCSymMatrix>;
00380 QCMANAGER_ITER_METH_EXPL_INST(void QCHamiltonGtr<QCSymMatrix>::computeElems);
00381
00382 QCPARAMETER_METH_EXPL_INST_PARAM(void QCHamiltonGtr<QCSymMatrix>::completeElems,
00383 SEVEN_PARAMS(const QCSubDomain& QCRestrict,
00384 const QCSubDomain& QCRestrict,
00385 QCRepInterElec& QCRestrict,
00386 QCSymMatrix& QCRestrict,
00387 const QCGeneralData&,
00388 int,
00389 int));