修订版 | 4932e7b9f31bf2e15fc93a21b24ac9e2d8d95e68 (tree) |
---|---|
时间 | 2013-05-23 10:32:52 |
作者 | Mikiya Fujii <mikiya.fujii@gmai...> |
Commiter | Mikiya Fujii |
Fast algorithm of ZindoS::GetAuxiliaryKNRKRElement is compelted. #31221
git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1350 1136aad2-a195-0410-b898-f5ea1d11b9d8
@@ -3344,7 +3344,6 @@ double ZindoS::GetKRElement(int moI, int moJ, int moK, int moL) const{ | ||
3344 | 3344 | double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{ |
3345 | 3345 | double value = 0.0; |
3346 | 3346 | |
3347 | -/* | |
3348 | 3347 | // Fast algorith, but this is not easy to read. |
3349 | 3348 | // Slow algorithm is alos written below. |
3350 | 3349 | for(int A=0; A<this->molecule->GetNumberAtoms(); A++){ |
@@ -3352,196 +3351,84 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons | ||
3352 | 3351 | int firstAOIndexA = atomA.GetFirstAOIndex(); |
3353 | 3352 | int lastAOIndexA = atomA.GetLastAOIndex(); |
3354 | 3353 | |
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 | - OrbitalType orbitalMu = atomA.GetValence(muOffSet); | |
3366 | - int lambda = mu; | |
3367 | - for(int nu=mu+1; nu<=lastAOIndexA; nu++){ | |
3368 | - int nuOffSet = nu - firstAOIndexA; | |
3369 | - OrbitalType orbitalNu = atomA.GetValence(nuOffSet); | |
3370 | - int sigma = nu; | |
3371 | - double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0, | |
3372 | - tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0, | |
3373 | - tmpMN07 = 0.0, tmpMN08 = 0.0, tmpMN09 = 0.0, | |
3374 | - tmpMN10 = 0.0, tmpMN11 = 0.0, tmpMN12 = 0.0, | |
3375 | - tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0, | |
3376 | - tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0, | |
3377 | - tmpMN19 = 0.0, tmpMN20 = 0.0, tmpMN21 = 0.0, | |
3378 | - tmpMN22 = 0.0, tmpMN23 = 0.0, tmpMN24 = 0.0; | |
3379 | - tmpMN01 = 4.0 | |
3380 | - *this->fockMatrix[moI][mu] | |
3381 | - *this->fockMatrix[moJ][nu]; | |
3382 | - tmpMN02 = 4.0 | |
3383 | - *this->fockMatrix[moK][mu] | |
3384 | - *this->fockMatrix[moL][nu]; | |
3385 | - tmpMN03 = this->fockMatrix[moI][mu] | |
3386 | - *this->fockMatrix[moK][nu]; | |
3387 | - tmpMN04 = this->fockMatrix[moJ][mu] | |
3388 | - *this->fockMatrix[moL][nu]; | |
3389 | - tmpMN05 = this->fockMatrix[moI][mu] | |
3390 | - *this->fockMatrix[moL][nu]; | |
3391 | - tmpMN06 = this->fockMatrix[moJ][mu] | |
3392 | - *this->fockMatrix[moK][nu]; | |
3393 | - tmpMN07 = 4.0 | |
3394 | - *this->fockMatrix[moI][mu] | |
3395 | - *this->fockMatrix[moJ][nu]; | |
3396 | - tmpMN08 = 4.0 | |
3397 | - *this->fockMatrix[moK][mu] | |
3398 | - *this->fockMatrix[moL][nu]; | |
3399 | - tmpMN09 = this->fockMatrix[moI][mu] | |
3400 | - *this->fockMatrix[moK][nu]; | |
3401 | - tmpMN10 = this->fockMatrix[moJ][mu] | |
3402 | - *this->fockMatrix[moL][nu]; | |
3403 | - tmpMN11 = this->fockMatrix[moI][mu] | |
3404 | - *this->fockMatrix[moL][nu]; | |
3405 | - tmpMN12 = this->fockMatrix[moJ][mu] | |
3406 | - *this->fockMatrix[moK][nu]; | |
3407 | - tmpMN13 = 4.0 | |
3408 | - *this->fockMatrix[moI][nu] | |
3354 | + for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){ | |
3355 | + OrbitalType orbitalMu = atomA.GetValence(mu - firstAOIndexA); | |
3356 | + double tmpM01 = this->fockMatrix[moI][mu] | |
3409 | 3357 | *this->fockMatrix[moJ][mu]; |
3410 | - tmpMN14 = 4.0 | |
3411 | - *this->fockMatrix[moK][nu] | |
3358 | + double tmpM02 = this->fockMatrix[moI][mu] | |
3412 | 3359 | *this->fockMatrix[moL][mu]; |
3413 | - tmpMN15 = this->fockMatrix[moI][nu] | |
3360 | + double tmpM03 = this->fockMatrix[moJ][mu] | |
3414 | 3361 | *this->fockMatrix[moK][mu]; |
3415 | - tmpMN16 = this->fockMatrix[moJ][nu] | |
3416 | - *this->fockMatrix[moL][mu]; | |
3417 | - tmpMN17 = this->fockMatrix[moI][nu] | |
3418 | - *this->fockMatrix[moL][mu]; | |
3419 | - tmpMN18 = this->fockMatrix[moJ][nu] | |
3362 | + double tmpM04 = this->fockMatrix[moI][mu] | |
3420 | 3363 | *this->fockMatrix[moK][mu]; |
3421 | - tmpMN19 = 4.0 | |
3422 | - *this->fockMatrix[moI][nu] | |
3423 | - *this->fockMatrix[moJ][mu]; | |
3424 | - tmpMN20 = 4.0 | |
3425 | - *this->fockMatrix[moK][nu] | |
3426 | - *this->fockMatrix[moL][mu]; | |
3427 | - tmpMN21 = this->fockMatrix[moI][nu] | |
3428 | - *this->fockMatrix[moK][mu]; | |
3429 | - tmpMN22 = this->fockMatrix[moJ][nu] | |
3364 | + double tmpM05 = this->fockMatrix[moJ][mu] | |
3430 | 3365 | *this->fockMatrix[moL][mu]; |
3431 | - tmpMN23 = this->fockMatrix[moI][nu] | |
3366 | + double tmpM06 = this->fockMatrix[moK][mu] | |
3432 | 3367 | *this->fockMatrix[moL][mu]; |
3433 | - tmpMN24 = this->fockMatrix[moJ][nu] | |
3434 | - *this->fockMatrix[moK][mu]; | |
3435 | 3368 | |
3436 | - { | |
3437 | - double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0, | |
3438 | - tmpMNL05 = 0.0, tmpMNL06 = 0.0, tmpMNL07 = 0.0, tmpMNL08 = 0.0, | |
3439 | - tmpMNL09 = 0.0, tmpMNL10 = 0.0, tmpMNL11 = 0.0, tmpMNL12 = 0.0, | |
3440 | - tmpMNL13 = 0.0, tmpMNL14 = 0.0, tmpMNL15 = 0.0, tmpMNL16 = 0.0, | |
3441 | - tmpMNL17 = 0.0, tmpMNL18 = 0.0, tmpMNL19 = 0.0, tmpMNL20 = 0.0, | |
3442 | - tmpMNL21 = 0.0, tmpMNL22 = 0.0, tmpMNL23 = 0.0, tmpMNL24 = 0.0; | |
3443 | - tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda]*this->fockMatrix[moL][sigma]; | |
3444 | - tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda]*this->fockMatrix[moJ][sigma]; | |
3445 | - tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda]*this->fockMatrix[moL][sigma]; | |
3446 | - tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda]*this->fockMatrix[moK][sigma]; | |
3447 | - tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda]*this->fockMatrix[moK][sigma]; | |
3448 | - tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda]*this->fockMatrix[moL][sigma]; | |
3449 | - tmpMNL07 = tmpMN07*this->fockMatrix[moL][lambda]*this->fockMatrix[moK][sigma]; | |
3450 | - tmpMNL08 = tmpMN08*this->fockMatrix[moJ][lambda]*this->fockMatrix[moI][sigma]; | |
3451 | - tmpMNL09 = tmpMN09*this->fockMatrix[moL][lambda]*this->fockMatrix[moJ][sigma]; | |
3452 | - tmpMNL10 = tmpMN10*this->fockMatrix[moK][lambda]*this->fockMatrix[moI][sigma]; | |
3453 | - tmpMNL11 = tmpMN11*this->fockMatrix[moK][lambda]*this->fockMatrix[moJ][sigma]; | |
3454 | - tmpMNL12 = tmpMN12*this->fockMatrix[moL][lambda]*this->fockMatrix[moI][sigma]; | |
3455 | - tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda]*this->fockMatrix[moL][sigma]; | |
3456 | - tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda]*this->fockMatrix[moJ][sigma]; | |
3457 | - tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda]*this->fockMatrix[moL][sigma]; | |
3458 | - tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda]*this->fockMatrix[moK][sigma]; | |
3459 | - tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda]*this->fockMatrix[moK][sigma]; | |
3460 | - tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda]*this->fockMatrix[moL][sigma]; | |
3461 | - tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda]*this->fockMatrix[moK][sigma]; | |
3462 | - tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda]*this->fockMatrix[moI][sigma]; | |
3463 | - tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda]*this->fockMatrix[moJ][sigma]; | |
3464 | - tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda]*this->fockMatrix[moI][sigma]; | |
3465 | - tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda]*this->fockMatrix[moJ][sigma]; | |
3466 | - tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda]*this->fockMatrix[moI][sigma]; | |
3467 | - | |
3468 | - tmpMNL01 +=-tmpMNL03 - tmpMNL06 + tmpMNL13 - tmpMNL15 - tmpMNL18; | |
3469 | - tmpMNL02 += tmpMNL14; | |
3470 | - tmpMNL04 += tmpMNL05 + tmpMNL16 + tmpMNL17; | |
3471 | - tmpMNL07 += tmpMNL19; | |
3472 | - tmpMNL08 +=-tmpMNL10 - tmpMNL12 + tmpMNL20 - tmpMNL22 - tmpMNL24; | |
3473 | - tmpMNL09 += tmpMNL11 + tmpMNL21 + tmpMNL23; | |
3474 | - { | |
3475 | - double tmpValue = 0.0; | |
3476 | - tmpValue += tmpMNL01; | |
3477 | - tmpValue += tmpMNL02; | |
3478 | - tmpValue -= tmpMNL04; | |
3479 | - tmpValue += tmpMNL07; | |
3480 | - tmpValue += tmpMNL08; | |
3481 | - tmpValue -= tmpMNL09; | |
3482 | - | |
3483 | - double gamma = 0.5*this->GetExchangeInt(orbitalMu, orbitalNu, atomA); | |
3484 | - value += tmpValue*gamma; | |
3485 | - } | |
3486 | - } | |
3487 | - } | |
3488 | - } | |
3489 | - } | |
3490 | - | |
3491 | - // (mu==nu && lambda==sigma) && (A==B || A!=B) | |
3492 | - for(int B=A; B<this->molecule->GetNumberAtoms(); B++){ | |
3493 | - const Atom& atomB = *this->molecule->GetAtom(B); | |
3494 | - int firstAOIndexB = atomB.GetFirstAOIndex(); | |
3495 | - int lastAOIndexB = atomB.GetLastAOIndex(); | |
3369 | + // A=B && (mu==lambda && nu==sigma for (mu nu|lamba sigma)) | |
3370 | + for(int nu=mu+1; nu<=lastAOIndexA; nu++){ | |
3371 | + OrbitalType orbitalNu = atomA.GetValence(nu - firstAOIndexA); | |
3372 | + double tmpValue = 0.0; | |
3373 | + tmpValue -= 4.0*tmpM01 | |
3374 | + *this->fockMatrix[moK][nu] | |
3375 | + *this->fockMatrix[moL][nu]; | |
3376 | + tmpValue += 6.0*tmpM02 | |
3377 | + *this->fockMatrix[moJ][nu] | |
3378 | + *this->fockMatrix[moK][nu]; | |
3379 | + tmpValue += 6.0*tmpM03 | |
3380 | + *this->fockMatrix[moL][nu] | |
3381 | + *this->fockMatrix[moI][nu]; | |
3382 | + tmpValue += 6.0*tmpM04 | |
3383 | + *this->fockMatrix[moJ][nu] | |
3384 | + *this->fockMatrix[moL][nu]; | |
3385 | + tmpValue += 6.0*tmpM05 | |
3386 | + *this->fockMatrix[moI][nu] | |
3387 | + *this->fockMatrix[moK][nu]; | |
3388 | + tmpValue -= 4.0*tmpM06 | |
3389 | + *this->fockMatrix[moI][nu] | |
3390 | + *this->fockMatrix[moJ][nu]; | |
3391 | + | |
3392 | + double gamma = 0.5*this->GetExchangeInt(orbitalMu, orbitalNu, atomA); | |
3393 | + value += tmpValue*gamma; | |
3394 | + } | |
3395 | + | |
3396 | + // (A==B || A!=B) && (mu==nu && lambda==sigma for (mu nu|lamba sigma)) | |
3397 | + for(int B=A; B<this->molecule->GetNumberAtoms(); B++){ | |
3398 | + const Atom& atomB = *this->molecule->GetAtom(B); | |
3399 | + int firstAOIndexB = atomB.GetFirstAOIndex(); | |
3400 | + int lastAOIndexB = atomB.GetLastAOIndex(); | |
3496 | 3401 | |
3497 | - for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){ | |
3498 | - int muOffSet = mu - firstAOIndexA; | |
3499 | - double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0, | |
3500 | - tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0; | |
3501 | - tmpMN01 = 4.0 | |
3502 | - *this->fockMatrix[moI][mu] | |
3503 | - *this->fockMatrix[moJ][mu]; | |
3504 | - tmpMN02 = 4.0 | |
3505 | - *this->fockMatrix[moK][mu] | |
3506 | - *this->fockMatrix[moL][mu]; | |
3507 | - tmpMN03 = this->fockMatrix[moI][mu] | |
3508 | - *this->fockMatrix[moK][mu]; | |
3509 | - tmpMN04 = this->fockMatrix[moJ][mu] | |
3510 | - *this->fockMatrix[moL][mu]; | |
3511 | - tmpMN05 = this->fockMatrix[moI][mu] | |
3512 | - *this->fockMatrix[moL][mu]; | |
3513 | - tmpMN06 = this->fockMatrix[moJ][mu] | |
3514 | - *this->fockMatrix[moK][mu]; | |
3515 | 3402 | for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){ |
3516 | - int lambdaOffSet = lambda - firstAOIndexB; | |
3517 | - double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0, | |
3518 | - tmpMNL05 = 0.0, tmpMNL06 = 0.0; | |
3519 | - tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda]; | |
3520 | - tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda]; | |
3521 | - tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda]; | |
3522 | - tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda]; | |
3523 | - tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda]; | |
3524 | - tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda]; | |
3525 | - tmpMNL01 -= tmpMNL03 + tmpMNL06; | |
3526 | - tmpMNL04 += tmpMNL05; | |
3403 | + OrbitalType orbitalLambda = atomB.GetValence(lambda - firstAOIndexB); | |
3527 | 3404 | double tmpValue = 0.0; |
3528 | - tmpValue += tmpMNL01*this->fockMatrix[moL][lambda]; | |
3529 | - tmpValue += tmpMNL02*this->fockMatrix[moJ][lambda]; | |
3530 | - tmpValue -= tmpMNL04*this->fockMatrix[moK][lambda]; | |
3405 | + tmpValue += 4.0*tmpM01 | |
3406 | + *this->fockMatrix[moK][lambda] | |
3407 | + *this->fockMatrix[moL][lambda]; | |
3408 | + tmpValue -= tmpM02 | |
3409 | + *this->fockMatrix[moJ][lambda] | |
3410 | + *this->fockMatrix[moK][lambda]; | |
3411 | + tmpValue -= tmpM03 | |
3412 | + *this->fockMatrix[moI][lambda] | |
3413 | + *this->fockMatrix[moL][lambda]; | |
3414 | + tmpValue -= tmpM04 | |
3415 | + *this->fockMatrix[moJ][lambda] | |
3416 | + *this->fockMatrix[moL][lambda]; | |
3417 | + tmpValue -= tmpM05 | |
3418 | + *this->fockMatrix[moI][lambda] | |
3419 | + *this->fockMatrix[moK][lambda]; | |
3420 | + tmpValue += 4.0*tmpM06 | |
3421 | + *this->fockMatrix[moI][lambda] | |
3422 | + *this->fockMatrix[moJ][lambda]; | |
3531 | 3423 | double gamma = 0.0; |
3532 | 3424 | if(A!=B){ |
3533 | - const double rAB = this->molecule->GetDistanceAtoms(atomA, atomB); | |
3534 | - OrbitalType orbitalMu = atomA.GetValence(muOffSet); | |
3535 | - OrbitalType orbitalLambda = atomB.GetValence(lambdaOffSet); | |
3536 | 3425 | gamma = this->GetNishimotoMatagaTwoEleInt(atomA, |
3537 | 3426 | orbitalMu, |
3538 | 3427 | atomB, |
3539 | 3428 | orbitalLambda, |
3540 | - rAB); | |
3429 | + this->molecule->GetDistanceAtoms(atomA, atomB)); | |
3541 | 3430 | } |
3542 | 3431 | else{ |
3543 | - OrbitalType orbitalMu = atomA.GetValence(muOffSet); | |
3544 | - OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet); | |
3545 | 3432 | gamma = 0.5*this->GetCoulombInt(orbitalMu, orbitalLambda, atomA); |
3546 | 3433 | } |
3547 | 3434 | value += tmpValue*gamma; |
@@ -3550,8 +3437,8 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons | ||
3550 | 3437 | } |
3551 | 3438 | } |
3552 | 3439 | // End of the fast algorith. |
3553 | -*/ | |
3554 | 3440 | |
3441 | +/* | |
3555 | 3442 | // slow algorithm |
3556 | 3443 | value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL, |
3557 | 3444 | *this->molecule, |
@@ -3562,7 +3449,7 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons | ||
3562 | 3449 | -1.0*this->GetMolecularIntegralElement(moI, moL, moJ, moK, |
3563 | 3450 | *this->molecule, |
3564 | 3451 | this->fockMatrix, NULL); |
3565 | - | |
3452 | +*/ | |
3566 | 3453 | return value; |
3567 | 3454 | } |
3568 | 3455 |