00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "QCCommon.hpp"
00021 #include "QCOverlapGtr.hpp"
00022 #include "QCSymMatrix.hpp"
00023 #include "QCManager.hpp"
00024 #include "QCMndo.hpp"
00025 #include "QCMndoParam.hpp"
00026 #include "QCAm1.hpp"
00027 #include "QCAm1Param.hpp"
00028 #include "QCPm3.hpp"
00029 #include "QCPm3Param.hpp"
00030 #include "QCGlobalSystem.hpp"
00031 #include "QCMDSystem.hpp"
00032 #include "QCDistMDSystem.hpp"
00033 #include "QCSCFAlgo.hpp"
00034 #include "QCDCAlgo.hpp"
00035
00036
00037
00038
00039
00040
00041
00042
00046 template <class TPMatrix>
00047 QCFloat *QCOverlapGtr<TPMatrix>::ExpHalfVector = NULL;
00048
00049
00050
00054 template <class TPMatrix>
00055 QCOverlapGtr<TPMatrix>::QCOverlapGtr (void) :
00056 QCMatElemGtr<TPMatrix> ()
00057 {}
00058
00059
00060
00061
00065 template <class TPMatrix>
00066 QCOverlapGtr<TPMatrix>::QCOverlapGtr (const int& dim) :
00067 QCMatElemGtr<TPMatrix> (dim)
00068 {}
00069
00070
00071
00072
00076 template <class TPMatrix>
00077 QCOverlapGtr<TPMatrix>::~QCOverlapGtr (void)
00078 {}
00079
00080
00081
00082
00086 template <class TPMatrix>
00087 template <class TPManager>
00088 void
00089 QCOverlapGtr<TPMatrix>::FillExpHalfVector (TPManager& manager)
00090 {
00091
00092 static int* quantNumNVector;
00093 static QCFloat* dzetaVector;
00094
00095 typedef typename TPManager::TModel TModel;
00096
00097 const typename TModel::TParam *QCRestrict params = manager.getParameters();
00098
00099
00100 quantNumNVector = new int [TModel::NbMaxAtomTypes];
00101 dzetaVector = new QCFloat [TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes];
00102 ExpHalfVector = new QCFloat [TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes];
00103
00104
00105 memset(quantNumNVector, 0,
00106 TModel::NbMaxAtomTypes * sizeof(int) );
00107 memset(dzetaVector, 0,
00108 TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes * sizeof(QCFloat) );
00109 memset(ExpHalfVector, 0,
00110 TModel::NbMaxOrbTypes * TModel::NbMaxAtomTypes * sizeof(QCFloat) );
00111
00112
00113
00114 FillQuantNumNVector(params, manager.getNbAtomTypes(), quantNumNVector);
00115 FillDzetaVector (params, manager.getNbAtomTypes(),
00116 TModel::NbMaxOrbTypes, dzetaVector);
00117
00123 QCFortran(fillexhalf)
00124 (quantNumNVector, dzetaVector, ExpHalfVector);
00125
00126 if (quantNumNVector) { delete [] quantNumNVector; }
00127 if (dzetaVector) { delete [] dzetaVector; }
00128 }
00129
00130
00131
00132
00136 template <class TPMatrix>
00137 void
00138 QCOverlapGtr<TPMatrix>::DestroyExpHalfVector (void) {
00139 if (ExpHalfVector)
00140 {
00141 delete [] ExpHalfVector;
00142 ExpHalfVector = NULL;
00143 }
00144 }
00145
00146
00147
00148
00152 template <class TPMatrix>
00153 template <class TPManager>
00154 void
00155 QCOverlapGtr<TPMatrix>::computeElems (TPManager& manager,
00156 typename TPManager::TSystem::QCIterator& workingSystem,
00157
00158
00159
00160 const QCPoint3D& shift)
00161 {
00162
00163 typedef typename TPManager::TModel TModel;
00164 typedef typename TPManager::TModel::TParam TParam;
00165
00170 TModel& QCRestrict model = manager.getModel();
00171 const TParam * QCRestrict params = manager.getParameters();
00172
00177 QCMatrix& QCRestrict overlapSAB = model.getSpdWorkingAB();
00178
00179
00180 computeElems(params,
00181 workingSystem->getSystem(),
00182 overlapSAB,
00183 TModel::NbMaxOrbTypes,
00184 shift);
00185 }
00186
00187
00188
00189
00190
00194 template <class TPMatrix>
00195 template <class TPParam>
00196 void
00197 QCOverlapGtr<TPMatrix>::computeElems (const TPParam * QCRestrict parameters,
00198 const QCSystem& QCRestrict system,
00199 QCMatrix& QCRestrict overlapSAB,
00200 int nbMaxOrbTypes,
00201 const QCPoint3D& QCRestrict shift) {
00202
00207 int firstAOofA, firstAOofB;
00208
00212 const TPParam *QCRestrict paramA;
00213 const TPParam *QCRestrict paramB;
00214
00215
00216 QCPoint3D coordsA, coordsB;
00217
00218
00223 for (int i = 1; i < system.getNbAtoms(); ++i) {
00224
00225 system.getPointAt(i, coordsA);
00226 paramA = ¶meters[system.getParamIndexAt(i)];
00227 firstAOofA = system.getFirstAOAt(i);
00228
00229 for (int mu = firstAOofA;
00230 mu < (firstAOofA + paramA->getNbAO());
00231 mu++)
00232 {
00233
00234 this->matrix[mu][mu] = QC_ONE;
00235 }
00236
00237 for (int j = 0; j < i; j++) {
00238
00239 system.getPointAt(j, coordsB);
00240 paramB = ¶meters[system.getParamIndexAt(j)];
00241 firstAOofB = system.getFirstAOAt(j);
00242
00243 computeElemsAux ( coordsA,
00244 *paramA,
00245 (paramA->getAtomicNbZ()-1) * nbMaxOrbTypes,
00246 coordsB,
00247 *paramB,
00248 (paramB->getAtomicNbZ()-1) * nbMaxOrbTypes,
00249 overlapSAB,
00250 shift );
00251
00258 this->matrix.insertSubMatrix( overlapSAB,
00259 firstAOofA, paramA->getNbAO(),
00260 firstAOofB, paramB->getNbAO() );
00261
00262 }
00263 }
00264 #ifdef DEBUG
00265
00266 #endif
00267 }
00268
00269
00270
00271
00272
00276 template <class TPMatrix>
00277 template <class TPParam>
00278 void
00279 QCOverlapGtr<TPMatrix>::computeElemsAux (const QCPoint3D& QCRestrict coordsA,
00280 const TPParam& QCRestrict paramA,
00281 int expHalfIndexA,
00282 const QCPoint3D& QCRestrict coordsB,
00283 const TPParam& QCRestrict paramB,
00284 int expHalfIndexB,
00285 QCMatrix& QCRestrict overlapSAB,
00286 const QCPoint3D& shift)
00287 {
00288
00289
00290 QCPoint3D coordsdB;
00291
00292 coordsdB[COORDX] = coordsB[COORDX] + shift[COORDX];
00293 coordsdB[COORDY] = coordsB[COORDY] + shift[COORDY];
00294 coordsdB[COORDZ] = coordsB[COORDZ] + shift[COORDZ];
00295
00296 QCFloat rAB = sqrt( QCPow<2>(coordsA[COORDX] - coordsdB[COORDX]) +
00297 QCPow<2>(coordsA[COORDY] - coordsdB[COORDY]) +
00298 QCPow<2>(coordsA[COORDZ] - coordsdB[COORDZ]) );
00299
00300 QCFortran(overlp)
00301 (
00302
00303
00304 paramB.getAtomicNbZ(),
00305 paramB.getNbAO(),
00306 paramB.getQuantumNumberN(),
00307 paramB.getDzetaS(),
00308 paramB.getDzetaP(),
00309 paramB.getDzetaD(),
00310 coordsdB[COORDX],
00311 coordsdB[COORDY],
00312 coordsdB[COORDZ],
00313 &ExpHalfVector[expHalfIndexB],
00314 paramA.getAtomicNbZ(),
00315 paramA.getNbAO(),
00316 paramA.getQuantumNumberN(),
00317 paramA.getDzetaS(),
00318 paramA.getDzetaP(),
00319 paramA.getDzetaD(),
00320 coordsA[COORDX],
00321 coordsA[COORDY],
00322 coordsA[COORDZ],
00323 &ExpHalfVector[expHalfIndexA],
00324 rAB,
00325 overlapSAB.getElems()
00326 );
00327
00328 }
00329
00330
00331
00332
00333
00337 template <class TPMatrix>
00338 template <class TPManager>
00339 void
00340 QCOverlapGtr<TPMatrix>::computeDerivatives (TPManager& manager,
00341 typename TPManager::TSystem::QCIterator& workingSystem,
00342 const QCPoint3D& shift,
00343 const QCFloat& QCRestrict inv2DeltaQ) {
00344
00345 QCOverlapGtr<QCSymMatrix>& QCRestrict
00346 workingIntg = *(workingSystem->getMatrices().getOverlapS());
00347
00348
00349 if (manager.getGeneralData().getDerivationType() == QC_NUMERICAL_DERIVATION) {
00350
00351
00352
00353
00354 QCPoint3D mshift;
00355 mshift -= shift;
00356 workingIntg.computeElems(manager, workingSystem, mshift);
00357
00358 #ifdef OUTPUT_DERIVATIVES2
00359 workingIntg.getMatrix().printInFile("DERIVATIVES_MS");
00360 #endif
00361
00362
00363
00364
00365 computeElems(manager, workingSystem, shift);
00366
00367 #ifdef OUTPUT_DERIVATIVES2
00368 this->matrix.printInFile("DERIVATIVES_S");
00369 #endif
00370
00371
00372 this->matrix.calculateNumericalDerivatives(workingIntg.getMatrix(),
00373 inv2DeltaQ);
00374
00375
00376 } else {
00377
00378 }
00379
00380 #ifdef OUTPUT_DERIVATIVES
00381 matrix.printInFile("DERIVATIVES");
00382 #endif
00383
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00401 template <class TPMatrix>
00402 template <class TPParam>
00403 void
00404 QCOverlapGtr<TPMatrix>::FillQuantNumNVector (const TPParam *QCRestrict params,
00405 const int& nbParameters,
00406 int* QCRestrict quantNumNVector)
00407 {
00408 for (int i=0; i < nbParameters; ++i) {
00409 quantNumNVector[params[i].getAtomicNbZ()-1] = params[i].getQuantumNumberN();
00410 }
00411 }
00412
00413
00414
00415
00421 template <class TPMatrix>
00422 template <class TPParam>
00423 void
00424 QCOverlapGtr<TPMatrix>::FillDzetaVector (const TPParam *QCRestrict params,
00425 const int& nbParameters,
00426 const int& nbMaxOrbTypes,
00427 QCFloat* QCRestrict dzetaVector)
00428 {
00429 int index;
00433 for (int i=0; i < nbParameters; ++i) {
00434
00435
00436 index = (params[i].getAtomicNbZ()-1) * nbMaxOrbTypes;
00437
00438
00439
00440 dzetaVector[index] = params[i].getDzetaS();
00441 dzetaVector[index+1] = params[i].getDzetaP();
00442
00443
00444 dzetaVector[index+2] = params[i].getDzetaD();
00445 }
00446 }
00447
00448
00449
00450
00454 template class QCOverlapGtr<QCSymMatrix>;
00455
00456 QCMANAGER_METH_EXPL_INST(void QCOverlapGtr<QCSymMatrix>::FillExpHalfVector);
00457
00458 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCOverlapGtr<QCSymMatrix>::computeElems,
00459 const QCPoint3D&);
00460
00461 QCMANAGER_ITER_METH_EXPL_INST_PARAM(void QCOverlapGtr<QCSymMatrix>::computeDerivatives,
00462 TWO_PARAMS(const QCPoint3D&,
00463 const QCFloat& QCRestrict));