00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "QCRepInterElec.hpp"
00021 #include "QCIntgReader.hpp"
00022 #include "QCIntgWriter.hpp"
00023 #include "QCFortran.hpp"
00024 #include "QCMndoParam.hpp"
00025 #include "QCAm1Param.hpp"
00026 #include "QCPm3Param.hpp"
00027 #include "QCSystem.hpp"
00028
00029
00030
00034 enum QCRepInterElecInd10 {
00035 s_s_10, s_px_10, s_py_10, s_pz_10,
00036 px_px_10, px_py_10, px_pz_10,
00037 py_py_10, py_pz_10,
00038 pz_pz_10
00039 };
00040 enum QCRepInterElecInd16 {
00041 s_s_16, s_px_16, s_py_16, s_pz_16,
00042 px_s_16, px_px_16, px_py_16, px_pz_16,
00043 py_s_16, py_px_16, py_py_16, py_pz_16,
00044 pz_s_16, pz_px_16, pz_py_16, pz_pz_16
00045 };
00046
00047
00048
00049
00053 QCRepInterElec::QCRepInterElec (const int& dim) :
00054 QCMatElemGtr<QCMatrix>(dim),
00055 intgWriter(NULL),
00056 intgReader(NULL),
00057 intgTmpFile("")
00058 {}
00059
00060
00061
00065 QCRepInterElec::~QCRepInterElec (void) {
00066 if (intgWriter) { delete intgWriter; intgWriter = NULL; }
00067 if (intgReader) { delete intgReader; intgReader = NULL; }
00068
00069 if (intgTmpFile != "") {
00070
00071
00072 remove(intgTmpFile.c_str());
00073 }
00074 }
00075
00076
00077
00081 QCIntgWriter *
00082 QCRepInterElec::allocIntgWriter (void) {
00083 if (!intgWriter) {
00084 intgWriter = new QCIntgWriter;
00085 }
00086 return intgWriter;
00087 }
00088
00089
00090
00094 QCIntgReader *
00095 QCRepInterElec::allocIntgReader (void) {
00096 if (!intgReader) {
00097 intgReader = new QCIntgReader;
00098 }
00099 return intgReader;
00100 }
00101
00102
00103
00108 void
00109 QCRepInterElec::expandOverDim16 (const int& nbAOofA, const int& nbAOofB) {
00110
00111
00112 const int SSRepIndex = 0;
00113 const int compressedDim = 10;
00114
00115 if (nbAOofA == 1) {
00124 matrix[SSRepIndex][pz_pz_16] = matrix[SSRepIndex][pz_pz_10];
00125 matrix[SSRepIndex][pz_py_16] = matrix[SSRepIndex][py_pz_10];
00126 matrix[SSRepIndex][pz_px_16] = matrix[SSRepIndex][px_pz_10];
00127 matrix[SSRepIndex][pz_s_16] = matrix[SSRepIndex][s_pz_10];
00133 matrix[SSRepIndex][py_py_16] = matrix[SSRepIndex][py_py_10];
00134 matrix[SSRepIndex][px_py_16] = matrix[SSRepIndex][px_py_10];
00135 matrix[SSRepIndex][py_px_16] = matrix[SSRepIndex][px_py_10];
00136 matrix[SSRepIndex][px_px_16] = matrix[SSRepIndex][px_px_10];
00140 matrix[SSRepIndex][px_s_16] = matrix[SSRepIndex][s_px_16];
00141 matrix[SSRepIndex][py_s_16] = matrix[SSRepIndex][s_py_16];
00142 matrix[SSRepIndex][px_pz_16] = matrix[SSRepIndex][pz_px_16];
00143 matrix[SSRepIndex][py_pz_16] = matrix[SSRepIndex][pz_py_16];
00144
00145 } else if (nbAOofB == 1) {
00149 matrix[pz_pz_16][SSRepIndex] = matrix[pz_pz_10][SSRepIndex];
00150 matrix[pz_py_16][SSRepIndex] = matrix[py_pz_10][SSRepIndex];
00151 matrix[pz_px_16][SSRepIndex] = matrix[px_pz_10][SSRepIndex];
00152 matrix[pz_s_16][SSRepIndex] = matrix[s_pz_10][SSRepIndex];
00158 matrix[py_py_16][SSRepIndex] = matrix[py_py_10][SSRepIndex];
00159 matrix[px_py_16][SSRepIndex] = matrix[px_py_10][SSRepIndex];
00160 matrix[py_px_16][SSRepIndex] = matrix[px_py_10][SSRepIndex];
00161 matrix[px_px_16][SSRepIndex] = matrix[px_px_10][SSRepIndex];
00165 matrix[px_s_16][SSRepIndex] = matrix[s_px_16][SSRepIndex];
00166 matrix[py_s_16][SSRepIndex] = matrix[s_py_16][SSRepIndex];
00167 matrix[px_pz_16][SSRepIndex] = matrix[pz_px_16][SSRepIndex];
00168 matrix[py_pz_16][SSRepIndex] = matrix[pz_py_16][SSRepIndex];
00169
00170 } else {
00174 for (int i = 0; i < compressedDim; ++i) {
00181 matrix[i][pz_pz_16] = matrix[i][pz_pz_10];
00182 matrix[i][pz_py_16] = matrix[i][py_pz_10];
00183 matrix[i][pz_px_16] = matrix[i][px_pz_10];
00184 matrix[i][pz_s_16] = matrix[i][s_pz_10];
00190 matrix[i][py_py_16] = matrix[i][py_py_10];
00191 matrix[i][px_py_16] = matrix[i][px_py_10];
00192 matrix[i][py_px_16] = matrix[i][px_py_10];
00193 matrix[i][px_px_16] = matrix[i][px_px_10];
00197 matrix[i][px_s_16] = matrix[i][s_px_16];
00198 matrix[i][py_s_16] = matrix[i][s_py_16];
00199 matrix[i][px_pz_16] = matrix[i][pz_px_16];
00200 matrix[i][py_pz_16] = matrix[i][pz_py_16];
00201 }
00211 memcpy(matrix[pz_pz_16], matrix[pz_pz_10], matrix.getDim()*sizeof(QCFloat) );
00212 memcpy(matrix[pz_py_16], matrix[py_pz_10], matrix.getDim()*sizeof(QCFloat) );
00213 memcpy(matrix[pz_px_16], matrix[px_pz_10], matrix.getDim()*sizeof(QCFloat) );
00214 memcpy(matrix[pz_s_16], matrix[s_pz_10], matrix.getDim()*sizeof(QCFloat) );
00220 memcpy(matrix[py_py_16], matrix[py_py_10], matrix.getDim()*sizeof(QCFloat) );
00221 memcpy(matrix[px_py_16], matrix[px_py_10], matrix.getDim()*sizeof(QCFloat) );
00222 memcpy(matrix[py_px_16], matrix[px_py_10], matrix.getDim()*sizeof(QCFloat) );
00223 memcpy(matrix[px_px_16], matrix[px_px_10], matrix.getDim()*sizeof(QCFloat) );
00227 memcpy(matrix[px_s_16], matrix[s_px_16], matrix.getDim()*sizeof(QCFloat) );
00228 memcpy(matrix[py_s_16], matrix[s_py_16], matrix.getDim()*sizeof(QCFloat) );
00229 memcpy(matrix[px_pz_16], matrix[pz_px_16], matrix.getDim()*sizeof(QCFloat) );
00230 memcpy(matrix[py_pz_16], matrix[pz_py_16], matrix.getDim()*sizeof(QCFloat) );
00231 }
00232 }
00233
00234
00235
00236
00241 template <class TPParam>
00242 void
00243 QCRepInterElec::computeElems (const TPParam* QCRestrict parameterA,
00244 const TPParam* QCRestrict parameterB,
00245 const QCPoint3D& QCRestrict coordsA,
00246 const QCPoint3D& QCRestrict coordsB,
00247 QCFloat rAB,
00248 bool ssACalc,
00249 bool ssBCalc,
00250 bool compressedFormat,
00251 const QCIntgAcquisitionMethod& intgAcquisitionMethod,
00252 bool isItfirstCall)
00253 {
00254
00255 int nbComputedAOofA;
00256 int nbComputedAOofB;
00257 int compressedFormatInt;
00262 if (compressedFormat) {
00263 compressedFormatInt = 1;
00264
00265 } else {
00266 compressedFormatInt = 0;
00267 }
00268
00269 if (ssBCalc) {
00270 nbComputedAOofB = TPParam::getNbSOrb();
00271
00272 } else {
00273 nbComputedAOofB = parameterB->getNbAO();
00274 }
00275
00284 if (ssACalc) {
00285 nbComputedAOofA = TPParam::getNbSOrb();
00286
00287 } else {
00288 nbComputedAOofA = parameterA->getNbAO();
00289 }
00290
00297 if ( (intgAcquisitionMethod == QC_DIRECT_COMPUTATION) ||
00298 (intgAcquisitionMethod == QC_INDIRECT_STORAGE && isItfirstCall) ) {
00299
00300 QCFortran(repul)
00301 (compressedFormatInt,
00302
00303 parameterB->getAtomicNbZ(),
00304
00305
00306 nbComputedAOofB,
00307 parameterB->getRo0(),
00308 parameterB->getRo1(),
00309 parameterB->getRo2(),
00310 parameterB->getD1(),
00311 parameterB->getD2(),
00312 coordsB[COORDX],
00313 coordsB[COORDY],
00314 coordsB[COORDZ],
00315
00316 parameterA->getAtomicNbZ(),
00317 nbComputedAOofA,
00318 parameterA->getRo0(),
00319 parameterA->getRo1(),
00320 parameterA->getRo2(),
00321 parameterA->getD1(),
00322 parameterA->getD2(),
00323 coordsA[COORDX],
00324 coordsA[COORDY],
00325 coordsA[COORDZ],
00326 rAB,
00327 matrix.getElems() );
00328
00329 if (intgAcquisitionMethod == QC_INDIRECT_STORAGE) {
00330
00331
00332 if (nbComputedAOofA > nbComputedAOofB) {
00333 matrix.pack(QCPow<2>(nbComputedAOofA) );
00334 }
00335
00336
00337 intgWriter->storeIntegrals(matrix.getElems(),
00338 sizeof(QCFloat),
00339 QCPow<2>(nbComputedAOofA) * QCPow<2>(nbComputedAOofB) );
00340
00341 if (nbComputedAOofA > nbComputedAOofB) {
00342 matrix.unpack(QCPow<2>(nbComputedAOofA) );
00343 }
00344 }
00345
00346 } else {
00347
00348
00349
00350
00351 intgReader->loadIntegrals(matrix.getElems(),
00352 sizeof(QCFloat),
00353 QCPow<2>(nbComputedAOofA) * QCPow<2>(nbComputedAOofB) );
00354
00355
00356 if (nbComputedAOofA > nbComputedAOofB) {
00357 matrix.unpack(QCPow<2>(nbComputedAOofA) );
00358 }
00359 }
00360
00361
00362 }
00363
00364
00365
00366
00367
00371 template <class TPParam>
00372 void
00373 QCRepInterElec::computeDerivatives (const TPParam* QCRestrict parameterA,
00374 const TPParam* QCRestrict parameterB,
00375 QCRepInterElec& QCRestrict workingIntg,
00376 QCDerivationType derivationType,
00377 const QCPoint3D& QCRestrict coordsA,
00378 const QCPoint3D& QCRestrict coordsB,
00379 const QCPoint3D& shift,
00380 bool ssACalc,
00381 bool ssBCalc,
00382 const QCFloat& QCRestrict inv2DeltaQ
00383 #ifdef TRACE_COMPUTE_DERIVATIVES
00384 ,
00385 QCChrono * intgChr1Ptr = NULL,
00386 QCChrono * intgChr2Ptr = NULL,
00387 QCChrono * intgChr3Ptr = NULL
00388 #endif
00389 ) {
00390
00391 int nbAOofA, nbAOofB;
00392 QCPoint3D coordsBdelta;
00393 QCFloat rAB;
00394 const int SSRepIndex = 0;
00401 bool compressedFormat = true;
00402
00403 if (derivationType == QC_NUMERICAL_DERIVATION) {
00404
00405 nbAOofA = parameterA->getNbAO();
00406 nbAOofB = parameterB->getNbAO();
00407
00408 coordsBdelta = coordsB;
00409 coordsBdelta -= shift;
00410
00411 rAB = sqrt(QCPow<2>(coordsBdelta[COORDX] - coordsA[COORDX]) +
00412 QCPow<2>(coordsBdelta[COORDY] - coordsA[COORDY]) +
00413 QCPow<2>(coordsBdelta[COORDZ] - coordsA[COORDZ]) );
00414
00415 #ifdef TRACE_COMPUTE_DERIVATIVES
00416 if (intgChr1Ptr) { intgChr1Ptr->run(); }
00417 #endif // TRACE_COMPUTE_DERIVATIVES
00418
00419 workingIntg.computeElems(
00420
00421 parameterA,
00422 parameterB,
00423 coordsA,
00424 coordsBdelta,
00425 rAB,
00426
00427
00428
00429 ssACalc,
00430 ssBCalc,
00431 compressedFormat,
00432
00433
00434 QC_DIRECT_COMPUTATION);
00435
00436
00437 #ifdef TRACE_COMPUTE_DERIVATIVES
00438 if (intgChr1Ptr) { intgChr1Ptr->pause(); }
00439 #endif // TRACE_COMPUTE_DERIVATIVES
00440
00441 coordsBdelta = coordsB;
00442 coordsBdelta += shift;
00443
00444 rAB = sqrt(QCPow<2>(coordsBdelta[COORDX] - coordsA[COORDX]) +
00445 QCPow<2>(coordsBdelta[COORDY] - coordsA[COORDY]) +
00446 QCPow<2>(coordsBdelta[COORDZ] - coordsA[COORDZ]) );
00447
00448 #ifdef TRACE_COMPUTE_DERIVATIVES
00449 if (intgChr2Ptr) { intgChr2Ptr->run(); }
00450 #endif // TRACE_COMPUTE_DERIVATIVES
00451
00452 computeElems(
00453
00454 parameterA,
00455 parameterB,
00456 coordsA,
00457 coordsBdelta,
00458 rAB,
00459
00460
00461
00462 ssACalc,
00463 ssBCalc,
00464 compressedFormat,
00465
00466
00467 QC_DIRECT_COMPUTATION);
00468
00469 #ifdef TRACE_COMPUTE_DERIVATIVES
00470 if (intgChr2Ptr) { intgChr2Ptr->pause(); }
00471 if (intgChr3Ptr) { intgChr3Ptr->run(); }
00472 #endif // TRACE_COMPUTE_DERIVATIVES
00473
00474 if (ssACalc && ssBCalc) {
00475
00480 if (inv2DeltaQ != QC_ZERO) {
00481 matrix[SSRepIndex][SSRepIndex] =
00482 (matrix[SSRepIndex][SSRepIndex] - workingIntg.getMatrix()[SSRepIndex][SSRepIndex]) *
00483 inv2DeltaQ;
00484
00485 } else {
00486
00487 matrix[SSRepIndex][SSRepIndex] =
00488 (matrix[SSRepIndex][SSRepIndex] - workingIntg.getMatrix()[SSRepIndex][SSRepIndex]);
00489 }
00490
00491 } else if (ssACalc && !ssBCalc) {
00492
00493
00494 matrix.calculateNumericalDerivatives(workingIntg.getMatrix(),
00495 inv2DeltaQ,
00496 1,
00497 (nbAOofB * (nbAOofB+1) ) / 2);
00498 expandOverDim16(1, nbAOofB);
00499 }
00500
00501 #ifdef DEBUG_QC
00502 else if (!ssACalc && ssBCalc) {
00503
00504
00505 cout << "Pb in File " << __FILE__ << " at line " << __LINE__ << endl;
00506 abort();
00507 }
00508 #endif //DEBUG_QC
00509
00510 else {
00511
00512
00513 matrix.calculateNumericalDerivatives(workingIntg.getMatrix(),
00514 inv2DeltaQ,
00515 (nbAOofA * (nbAOofA+1) ) / 2 ,
00516 (nbAOofB * (nbAOofB+1) ) / 2);
00517 expandOverDim16(nbAOofA, nbAOofB);
00518 }
00519
00520 #ifdef TRACE_COMPUTE_DERIVATIVES
00521 if (intgChr3Ptr) { intgChr3Ptr->pause(); }
00522 #endif // TRACE_COMPUTE_DERIVATIVES
00523
00524
00525 } else {
00526
00527 }
00528
00529 }
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00544 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeElems,
00545 EIGHT_PARAMS(const QCPoint3D& QCRestrict,
00546 const QCPoint3D& QCRestrict,
00547 QCFloat,
00548 bool,
00549 bool,
00550 bool,
00551 const QCIntgAcquisitionMethod&,
00552 bool));
00553
00554 #ifndef TRACE_COMPUTE_DERIVATIVES
00555
00556 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeDerivatives,
00557 EIGHT_PARAMS(QCRepInterElec&,
00558 QCDerivationType,
00559 const QCPoint3D&,
00560 const QCPoint3D&,
00561 const QCPoint3D&,
00562 bool,
00563 bool,
00564 const QCFloat& QCRestrict));
00565
00566 #else // TRACE_COMPUTE_DERIVATIVES
00567 TWO_QCPARAMETERS_METH_EXPL_INST_PARAM(void QCRepInterElec::computeDerivatives,
00568 ELEVEN_PARAMS(QCRepInterElec&,
00569 QCDerivationType,
00570 const QCPoint3D&,
00571 const QCPoint3D,
00572 const QCPoint3D&,
00573 bool,
00574 bool,
00575 const QCFloat& QCRestrict
00576 QCChrono*,
00577 QCChrono*,
00578 QCChrono*));
00579
00580 #endif // TRACE_COMPUTE_DERIVATIVES