OSDN -- Desarrollar y descargar software de código abierto
Translation Status of Español translation status for es
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
No Tags

Frequently used words (click to add to your profile)

javac++androidlinuxc#windowsobjective-ccocoa誰得qtpythonphprubygameguibathyscaphec計画中(planning stage)翻訳omegatframeworktwitterdomtestvb.netdirectxゲームエンジンbtronarduinopreviewer

Commit MetaInfo

Revisión758257116ada46936474bc2d8291db3984121190 (tree)
Tiempo2013-05-22 21:03:27
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Fast algorithm of ZindoS::GetAuxiliaryKNRKRElement is tested. #31221

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1348 1136aad2-a195-0410-b898-f5ea1d11b9d8

Cambiar Resumen

Diferencia incremental

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -3352,29 +3352,36 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
33523352 int firstAOIndexA = atomA.GetFirstAOIndex();
33533353 int lastAOIndexA = atomA.GetLastAOIndex();
33543354
3355- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3356- int muOffSet = mu - firstAOIndexA;
3357- for(int nu=mu; nu<=lastAOIndexA; nu++){
3358- int nuOffSet = nu - firstAOIndexA;
3359- double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
3360- tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0,
3361- tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0,
3362- tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0;
3363- tmpMN01 = 4.0
3364- *this->fockMatrix[moI][mu]
3365- *this->fockMatrix[moJ][nu];
3366- tmpMN02 = 4.0
3367- *this->fockMatrix[moK][mu]
3368- *this->fockMatrix[moL][nu];
3369- tmpMN03 = this->fockMatrix[moI][mu]
3370- *this->fockMatrix[moK][nu];
3371- tmpMN04 = this->fockMatrix[moJ][mu]
3372- *this->fockMatrix[moL][nu];
3373- tmpMN05 = this->fockMatrix[moI][mu]
3374- *this->fockMatrix[moL][nu];
3375- tmpMN06 = this->fockMatrix[moJ][mu]
3376- *this->fockMatrix[moK][nu];
3377- if(mu != nu){
3355+ // A=B && (mu==lambda && nu==sigma)
3356+ //for(int B=A; B<A+1; B++){
3357+ {
3358+ int B=A;
3359+ const Atom& atomB = *this->molecule->GetAtom(B);
3360+ int firstAOIndexB = atomB.GetFirstAOIndex();
3361+ int lastAOIndexB = atomB.GetLastAOIndex();
3362+
3363+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3364+ int muOffSet = mu - firstAOIndexA;
3365+ for(int nu=mu+1; nu<=lastAOIndexA; nu++){
3366+ int nuOffSet = nu - firstAOIndexA;
3367+ double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
3368+ tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0,
3369+ tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0,
3370+ tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0;
3371+ tmpMN01 = 4.0
3372+ *this->fockMatrix[moI][mu]
3373+ *this->fockMatrix[moJ][nu];
3374+ tmpMN02 = 4.0
3375+ *this->fockMatrix[moK][mu]
3376+ *this->fockMatrix[moL][nu];
3377+ tmpMN03 = this->fockMatrix[moI][mu]
3378+ *this->fockMatrix[moK][nu];
3379+ tmpMN04 = this->fockMatrix[moJ][mu]
3380+ *this->fockMatrix[moL][nu];
3381+ tmpMN05 = this->fockMatrix[moI][mu]
3382+ *this->fockMatrix[moL][nu];
3383+ tmpMN06 = this->fockMatrix[moJ][mu]
3384+ *this->fockMatrix[moK][nu];
33783385 tmpMN13 = 4.0
33793386 *this->fockMatrix[moI][nu]
33803387 *this->fockMatrix[moJ][mu];
@@ -3389,14 +3396,9 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
33893396 *this->fockMatrix[moL][mu];
33903397 tmpMN18 = this->fockMatrix[moJ][nu]
33913398 *this->fockMatrix[moK][mu];
3392- }
33933399
3394- for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
3395- const Atom& atomB = *this->molecule->GetAtom(B);
3396- int firstAOIndexB = atomB.GetFirstAOIndex();
3397- int lastAOIndexB = atomB.GetLastAOIndex();
3398-
3399- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3400+ {
3401+ int lambda = mu;
34003402 int lambdaOffSet = lambda - firstAOIndexB;
34013403 double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0,
34023404 tmpMNL05 = 0.0, tmpMNL06 = 0.0, tmpMNL07 = 0.0, tmpMNL08 = 0.0,
@@ -3420,312 +3422,113 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
34203422 tmpMNL04 += tmpMNL05;
34213423 tmpMNL08 -= tmpMNL10 + tmpMNL12;
34223424 tmpMNL09 += tmpMNL11;
3423- if(mu != nu){
3424- tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda];
3425- tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda];
3426- tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda];
3427- tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda];
3428- tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda];
3429- tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda];
3430- tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda];
3431- tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda];
3432- tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda];
3433- tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda];
3434- tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda];
3435- tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda];
3436- tmpMNL13 -= tmpMNL15 + tmpMNL18;
3437- tmpMNL16 += tmpMNL17;
3438- tmpMNL20 -= tmpMNL22 + tmpMNL24;
3439- tmpMNL21 += tmpMNL23;
3440- tmpMNL01 += tmpMNL13;
3441- tmpMNL02 += tmpMNL14;
3442- tmpMNL04 += tmpMNL16;
3443- tmpMNL07 += tmpMNL19;
3444- tmpMNL08 += tmpMNL20;
3445- tmpMNL09 += tmpMNL21;
3446- }
3447- for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
3425+ tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda];
3426+ tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda];
3427+ tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda];
3428+ tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda];
3429+ tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda];
3430+ tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda];
3431+ tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda];
3432+ tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda];
3433+ tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda];
3434+ tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda];
3435+ tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda];
3436+ tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda];
3437+ tmpMNL13 -= tmpMNL15 + tmpMNL18;
3438+ tmpMNL16 += tmpMNL17;
3439+ tmpMNL20 -= tmpMNL22 + tmpMNL24;
3440+ tmpMNL21 += tmpMNL23;
3441+
3442+ tmpMNL01 += tmpMNL13;
3443+ tmpMNL02 += tmpMNL14;
3444+ tmpMNL04 += tmpMNL16;
3445+ tmpMNL07 += tmpMNL19;
3446+ tmpMNL08 += tmpMNL20;
3447+ tmpMNL09 += tmpMNL21;
3448+ {
3449+ int sigma = nu;
34483450 int sigmaOffSet = sigma - firstAOIndexB;
34493451 double tmpValue = 0.0;
34503452 tmpValue += tmpMNL01*this->fockMatrix[moL][sigma];
34513453 tmpValue += tmpMNL02*this->fockMatrix[moJ][sigma];
34523454 tmpValue -= tmpMNL04*this->fockMatrix[moK][sigma];
3453- if(lambda != sigma){
3454- tmpValue += tmpMNL07*this->fockMatrix[moK][sigma];
3455- tmpValue += tmpMNL08*this->fockMatrix[moI][sigma];
3456- tmpValue -= tmpMNL09*this->fockMatrix[moJ][sigma];
3457- }
3458- double gamma = 0.0;
3459- if(A!=B){
3460- gamma = this->twoElecTwoCore[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
3461- }
3462- else{
3463- if(mu==nu && lambda==sigma){
3464- OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3465- OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet);
3466- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3467- }
3468- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3469- OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3470- OrbitalType orbitalNu = atomA.GetValence(nuOffSet);
3471- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3472- }
3473- else{
3474- gamma = 0.0;
3475- }
3476- gamma *= 0.5;
3477- }
3455+ tmpValue += tmpMNL07*this->fockMatrix[moK][sigma];
3456+ tmpValue += tmpMNL08*this->fockMatrix[moI][sigma];
3457+ tmpValue -= tmpMNL09*this->fockMatrix[moJ][sigma];
3458+
3459+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3460+ OrbitalType orbitalNu = atomA.GetValence(nuOffSet);
3461+ double gamma = 0.5*this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
34783462 value += tmpValue*gamma;
34793463 }
34803464 }
34813465 }
34823466 }
34833467 }
3484- }
3485- // End of the fast algorith.
3486-*/
3487- /*
3488- // Algorithm using blas
3489- double** twoElec = NULL;
3490- double* twiceMoIJ = NULL;
3491- double* twiceMoIK = NULL;
3492- double* twiceMoIL = NULL;
3493- double* twiceMoKL = NULL;
3494- double* twiceMoJL = NULL;
3495- double* twiceMoJK = NULL;
3496- double* tmpVector = NULL;
3497- int numAOs = this->molecule->GetTotalNumberAOs();
3498- MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3499- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3500- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3501- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3502- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
3503- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
3504- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
3505- MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
3506- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3507- const Atom& atomA = *this->molecule->GetAtom(A);
3508- int firstAOIndexA = atomA.GetFirstAOIndex();
3509- int lastAOIndexA = atomA.GetLastAOIndex();
3510- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3511- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3512- twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
3513- twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
3514- twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
3515- }
3516- }
3517- }
3518-
3519- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3520- const Atom& atomB = *this->molecule->GetAtom(B);
3521- int firstAOIndexB = atomB.GetFirstAOIndex();
3522- int lastAOIndexB = atomB.GetLastAOIndex();
3523- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3524- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3525- twiceMoKL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
3526- twiceMoJL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
3527- twiceMoJK[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
3528- }
3529- }
3530- }
35313468
3532- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3533- const Atom& atomA = *this->molecule->GetAtom(A);
3534- int firstAOIndexA = atomA.GetFirstAOIndex();
3535- int lastAOIndexA = atomA.GetLastAOIndex();
3469+ // (mu==nu && lambda==sigma) && (A==B || A!=B)
35363470 for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
35373471 const Atom& atomB = *this->molecule->GetAtom(B);
35383472 int firstAOIndexB = atomB.GetFirstAOIndex();
35393473 int lastAOIndexB = atomB.GetLastAOIndex();
3540- double gamma = 0.0;
3541- if(A!=B){
3542- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3543- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3544- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3545- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3546- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3547- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3548- this->twoElecTwoCore[A]
3549- [B]
3550- [mu-firstAOIndexA]
3551- [nu-firstAOIndexA]
3552- [lambda-firstAOIndexB]
3553- [sigma-firstAOIndexB];
3554- }
3555- }
3556- }
3557- }
3558- }
3559- else{
3560- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3561- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3562- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3563- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3564- if(mu==nu && lambda==sigma){
3565- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3566- OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3567- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3568- }
3569- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3570- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3571- OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
3572- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3573- }
3574- else{
3575- gamma = 0.0;
3576- }
3577- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3578- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
3579- }
3580- }
3581- }
3582- }
3583- }
3584- }
3585- }
3586- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3587- twoElec,
3588- twiceMoKL,
3589- tmpVector);
3590- value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, tmpVector);
3591- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3592- twoElec,
3593- twiceMoJL,
3594- tmpVector);
3595- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, tmpVector);
3596- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3597- twoElec,
3598- twiceMoJK,
3599- tmpVector);
3600- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, tmpVector);
3601- MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3602- MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3603- MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3604- MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3605- MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
3606- MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
3607- MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
3608- MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
3609- // End of algorithm using blas
3610- */
36113474
3612- /*
3613- // Second algorithm using blas.
3614- // This algorithm uses DGEMM.
3615- double** twoElec = NULL;
3616- double* twiceMoIJ = NULL;
3617- double* twiceMoIK = NULL;
3618- double* twiceMoIL = NULL;
3619- double** twiceMoB = NULL;
3620- double** tmpMatrix = NULL;
3621- int numAOs = this->molecule->GetTotalNumberAOs();
3622- MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3623- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3624- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3625- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3626- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
3627- MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
3628- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3629- const Atom& atomA = *this->molecule->GetAtom(A);
3630- int firstAOIndexA = atomA.GetFirstAOIndex();
3631- int lastAOIndexA = atomA.GetLastAOIndex();
3632- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3633- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3634- twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
3635- twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
3636- twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
3637- }
3638- }
3639- }
3640-
3641- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3642- const Atom& atomB = *this->molecule->GetAtom(B);
3643- int firstAOIndexB = atomB.GetFirstAOIndex();
3644- int lastAOIndexB = atomB.GetLastAOIndex();
3645- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3646- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3647- twiceMoB[0][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
3648- twiceMoB[1][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
3649- twiceMoB[2][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
3650- }
3651- }
3652- }
3653-
3654- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3655- const Atom& atomA = *this->molecule->GetAtom(A);
3656- int firstAOIndexA = atomA.GetFirstAOIndex();
3657- int lastAOIndexA = atomA.GetLastAOIndex();
3658- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3659- const Atom& atomB = *this->molecule->GetAtom(B);
3660- int firstAOIndexB = atomB.GetFirstAOIndex();
3661- int lastAOIndexB = atomB.GetLastAOIndex();
3662- double gamma = 0.0;
3663- if(A!=B){
3664- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3665- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3666- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3667- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3668- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3669- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3670- this->twoElecTwoCore[A]
3671- [B]
3672- [mu-firstAOIndexA]
3673- [nu-firstAOIndexA]
3674- [lambda-firstAOIndexB]
3675- [sigma-firstAOIndexB];
3676- }
3677- }
3475+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3476+ int muOffSet = mu - firstAOIndexA;
3477+ double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
3478+ tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0;
3479+ tmpMN01 = 4.0
3480+ *this->fockMatrix[moI][mu]
3481+ *this->fockMatrix[moJ][mu];
3482+ tmpMN02 = 4.0
3483+ *this->fockMatrix[moK][mu]
3484+ *this->fockMatrix[moL][mu];
3485+ tmpMN03 = this->fockMatrix[moI][mu]
3486+ *this->fockMatrix[moK][mu];
3487+ tmpMN04 = this->fockMatrix[moJ][mu]
3488+ *this->fockMatrix[moL][mu];
3489+ tmpMN05 = this->fockMatrix[moI][mu]
3490+ *this->fockMatrix[moL][mu];
3491+ tmpMN06 = this->fockMatrix[moJ][mu]
3492+ *this->fockMatrix[moK][mu];
3493+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3494+ int lambdaOffSet = lambda - firstAOIndexB;
3495+ double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0,
3496+ tmpMNL05 = 0.0, tmpMNL06 = 0.0;
3497+ tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda];
3498+ tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda];
3499+ tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda];
3500+ tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda];
3501+ tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda];
3502+ tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda];
3503+ tmpMNL01 -= tmpMNL03 + tmpMNL06;
3504+ tmpMNL04 += tmpMNL05;
3505+ double tmpValue = 0.0;
3506+ tmpValue += tmpMNL01*this->fockMatrix[moL][lambda];
3507+ tmpValue += tmpMNL02*this->fockMatrix[moJ][lambda];
3508+ tmpValue -= tmpMNL04*this->fockMatrix[moK][lambda];
3509+ double gamma = 0.0;
3510+ if(A!=B){
3511+ const double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
3512+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3513+ OrbitalType orbitalLambda = atomB.GetValence(lambdaOffSet);
3514+ gamma = this->GetNishimotoMatagaTwoEleInt(atomA,
3515+ orbitalMu,
3516+ atomB,
3517+ orbitalLambda,
3518+ rAB);
36783519 }
3679- }
3680- }
3681- else{
3682- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3683- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3684- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3685- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3686- if(mu==nu && lambda==sigma){
3687- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3688- OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3689- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3690- }
3691- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3692- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3693- OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
3694- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3695- }
3696- else{
3697- gamma = 0.0;
3698- }
3699- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3700- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
3701- }
3702- }
3520+ else{
3521+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3522+ OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet);
3523+ gamma = 0.5*this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
37033524 }
3525+ value += tmpValue*gamma;
37043526 }
37053527 }
37063528 }
37073529 }
3708-
3709- MolDS_wrappers::Blas::GetInstance()->Dgemm(false, true, true,
3710- this->molecule->GetNumberAtoms()*dxy*dxy,
3711- 3,
3712- this->molecule->GetNumberAtoms()*dxy*dxy,
3713- 1.0,
3714- twoElec,
3715- twiceMoB,
3716- 0.0,
3717- tmpMatrix);
3718- value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
3719- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
3720- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
3721- MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3722- MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3723- MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3724- MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3725- MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
3726- MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
3727- // End of second algorithm using blas
3728- */
3530+ // End of the fast algorith.
3531+*/
37293532
37303533 // slow algorithm
37313534 value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL,