OSDN -- Open-Source-Software-Entwicklung und Downloads
Übersetzungsstatus von Deutsch translation status for de
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
Keine Tags

Frequently used words (click to add to your profile)

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

Commit MetaInfo

Revision8432e5e1c2771e7a80e141bedff28c173a672787 (tree)
Zeit2013-05-21 19:39:32
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

GetAuxiliaryKNRKRElement in ZindoS and Mndo is refactored. #31221

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

Ändern Zusammenfassung

Diff

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -2434,6 +2434,408 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs, dou
24342434 CartesianType_end);
24352435 }
24362436
2437+// see common term in eqs. (45) and (46) in [PT_1996],
2438+// that is, 4.0(ij|kl) - (ik|jl) - (il|jk).
2439+double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
2440+ double value = 0.0;
2441+
2442+ // Fast algorith, but this is not easy to read.
2443+ // Slow algorithm is alos written below.
2444+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
2445+ const Atom& atomA = *this->molecule->GetAtom(A);
2446+ int firstAOIndexA = atomA.GetFirstAOIndex();
2447+ int lastAOIndexA = atomA.GetLastAOIndex();
2448+
2449+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2450+ int muOffSet = mu - firstAOIndexA;
2451+ for(int nu=mu; nu<=lastAOIndexA; nu++){
2452+ int nuOffSet = nu - firstAOIndexA;
2453+ double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
2454+ tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0,
2455+ tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0,
2456+ tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0;
2457+ tmpMN01 = 4.0
2458+ *this->fockMatrix[moI][mu]
2459+ *this->fockMatrix[moJ][nu];
2460+ tmpMN02 = 4.0
2461+ *this->fockMatrix[moK][mu]
2462+ *this->fockMatrix[moL][nu];
2463+ tmpMN03 = this->fockMatrix[moI][mu]
2464+ *this->fockMatrix[moK][nu];
2465+ tmpMN04 = this->fockMatrix[moJ][mu]
2466+ *this->fockMatrix[moL][nu];
2467+ tmpMN05 = this->fockMatrix[moI][mu]
2468+ *this->fockMatrix[moL][nu];
2469+ tmpMN06 = this->fockMatrix[moJ][mu]
2470+ *this->fockMatrix[moK][nu];
2471+ if(mu != nu){
2472+ tmpMN13 = 4.0
2473+ *this->fockMatrix[moI][nu]
2474+ *this->fockMatrix[moJ][mu];
2475+ tmpMN14 = 4.0
2476+ *this->fockMatrix[moK][nu]
2477+ *this->fockMatrix[moL][mu];
2478+ tmpMN15 = this->fockMatrix[moI][nu]
2479+ *this->fockMatrix[moK][mu];
2480+ tmpMN16 = this->fockMatrix[moJ][nu]
2481+ *this->fockMatrix[moL][mu];
2482+ tmpMN17 = this->fockMatrix[moI][nu]
2483+ *this->fockMatrix[moL][mu];
2484+ tmpMN18 = this->fockMatrix[moJ][nu]
2485+ *this->fockMatrix[moK][mu];
2486+ }
2487+
2488+ for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
2489+ const Atom& atomB = *this->molecule->GetAtom(B);
2490+ int firstAOIndexB = atomB.GetFirstAOIndex();
2491+ int lastAOIndexB = atomB.GetLastAOIndex();
2492+
2493+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2494+ int lambdaOffSet = lambda - firstAOIndexB;
2495+ double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0,
2496+ tmpMNL05 = 0.0, tmpMNL06 = 0.0, tmpMNL07 = 0.0, tmpMNL08 = 0.0,
2497+ tmpMNL09 = 0.0, tmpMNL10 = 0.0, tmpMNL11 = 0.0, tmpMNL12 = 0.0,
2498+ tmpMNL13 = 0.0, tmpMNL14 = 0.0, tmpMNL15 = 0.0, tmpMNL16 = 0.0,
2499+ tmpMNL17 = 0.0, tmpMNL18 = 0.0, tmpMNL19 = 0.0, tmpMNL20 = 0.0,
2500+ tmpMNL21 = 0.0, tmpMNL22 = 0.0, tmpMNL23 = 0.0, tmpMNL24 = 0.0;
2501+ tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda];
2502+ tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda];
2503+ tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda];
2504+ tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda];
2505+ tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda];
2506+ tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda];
2507+ tmpMNL07 = tmpMN01*this->fockMatrix[moL][lambda];
2508+ tmpMNL08 = tmpMN02*this->fockMatrix[moJ][lambda];
2509+ tmpMNL09 = tmpMN03*this->fockMatrix[moL][lambda];
2510+ tmpMNL10 = tmpMN04*this->fockMatrix[moK][lambda];
2511+ tmpMNL11 = tmpMN05*this->fockMatrix[moK][lambda];
2512+ tmpMNL12 = tmpMN06*this->fockMatrix[moL][lambda];
2513+ tmpMNL01 -= tmpMNL03 + tmpMNL06;
2514+ tmpMNL04 += tmpMNL05;
2515+ tmpMNL08 -= tmpMNL10 + tmpMNL12;
2516+ tmpMNL09 += tmpMNL11;
2517+ if(mu != nu){
2518+ tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda];
2519+ tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda];
2520+ tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda];
2521+ tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda];
2522+ tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda];
2523+ tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda];
2524+ tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda];
2525+ tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda];
2526+ tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda];
2527+ tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda];
2528+ tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda];
2529+ tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda];
2530+ tmpMNL13 -= tmpMNL15 + tmpMNL18;
2531+ tmpMNL16 += tmpMNL17;
2532+ tmpMNL20 -= tmpMNL22 + tmpMNL24;
2533+ tmpMNL21 += tmpMNL23;
2534+ tmpMNL01 += tmpMNL13;
2535+ tmpMNL02 += tmpMNL14;
2536+ tmpMNL04 += tmpMNL16;
2537+ tmpMNL07 += tmpMNL19;
2538+ tmpMNL08 += tmpMNL20;
2539+ tmpMNL09 += tmpMNL21;
2540+ }
2541+ for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
2542+ int sigmaOffSet = sigma - firstAOIndexB;
2543+ double tmpValue = 0.0;
2544+ tmpValue += tmpMNL01*this->fockMatrix[moL][sigma];
2545+ tmpValue += tmpMNL02*this->fockMatrix[moJ][sigma];
2546+ tmpValue -= tmpMNL04*this->fockMatrix[moK][sigma];
2547+ if(lambda != sigma){
2548+ tmpValue += tmpMNL07*this->fockMatrix[moK][sigma];
2549+ tmpValue += tmpMNL08*this->fockMatrix[moI][sigma];
2550+ tmpValue -= tmpMNL09*this->fockMatrix[moJ][sigma];
2551+ }
2552+ double gamma = 0.0;
2553+ if(A!=B){
2554+ gamma = this->twoElecTwoCore[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
2555+ }
2556+ else{
2557+ if(mu==nu && lambda==sigma){
2558+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
2559+ OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet);
2560+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
2561+ }
2562+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
2563+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
2564+ OrbitalType orbitalNu = atomA.GetValence(nuOffSet);
2565+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
2566+ }
2567+ else{
2568+ gamma = 0.0;
2569+ }
2570+ gamma *= 0.5;
2571+ }
2572+ value += tmpValue*gamma;
2573+ }
2574+ }
2575+ }
2576+ }
2577+ }
2578+ }
2579+ // End of the fast algorith.
2580+
2581+ /*
2582+ // Algorithm using blas
2583+ double** twoElec = NULL;
2584+ double* twiceMoIJ = NULL;
2585+ double* twiceMoIK = NULL;
2586+ double* twiceMoIL = NULL;
2587+ double* twiceMoKL = NULL;
2588+ double* twiceMoJL = NULL;
2589+ double* twiceMoJK = NULL;
2590+ double* tmpVector = NULL;
2591+ int numAOs = this->molecule->GetTotalNumberAOs();
2592+ MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
2593+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
2594+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
2595+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
2596+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
2597+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
2598+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
2599+ MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
2600+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
2601+ const Atom& atomA = *this->molecule->GetAtom(A);
2602+ int firstAOIndexA = atomA.GetFirstAOIndex();
2603+ int lastAOIndexA = atomA.GetLastAOIndex();
2604+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2605+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2606+ twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
2607+ twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
2608+ twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
2609+ }
2610+ }
2611+ }
2612+
2613+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
2614+ const Atom& atomB = *this->molecule->GetAtom(B);
2615+ int firstAOIndexB = atomB.GetFirstAOIndex();
2616+ int lastAOIndexB = atomB.GetLastAOIndex();
2617+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2618+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2619+ twiceMoKL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
2620+ twiceMoJL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
2621+ twiceMoJK[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
2622+ }
2623+ }
2624+ }
2625+
2626+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
2627+ const Atom& atomA = *this->molecule->GetAtom(A);
2628+ int firstAOIndexA = atomA.GetFirstAOIndex();
2629+ int lastAOIndexA = atomA.GetLastAOIndex();
2630+ for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
2631+ const Atom& atomB = *this->molecule->GetAtom(B);
2632+ int firstAOIndexB = atomB.GetFirstAOIndex();
2633+ int lastAOIndexB = atomB.GetLastAOIndex();
2634+ double gamma = 0.0;
2635+ if(A!=B){
2636+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2637+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2638+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2639+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2640+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
2641+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
2642+ this->twoElecTwoCore[A]
2643+ [B]
2644+ [mu-firstAOIndexA]
2645+ [nu-firstAOIndexA]
2646+ [lambda-firstAOIndexB]
2647+ [sigma-firstAOIndexB];
2648+ }
2649+ }
2650+ }
2651+ }
2652+ }
2653+ else{
2654+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2655+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2656+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2657+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2658+ if(mu==nu && lambda==sigma){
2659+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2660+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
2661+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
2662+ }
2663+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
2664+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2665+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
2666+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
2667+ }
2668+ else{
2669+ gamma = 0.0;
2670+ }
2671+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
2672+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
2673+ }
2674+ }
2675+ }
2676+ }
2677+ }
2678+ }
2679+ }
2680+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
2681+ twoElec,
2682+ twiceMoKL,
2683+ tmpVector);
2684+ value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, tmpVector);
2685+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
2686+ twoElec,
2687+ twiceMoJL,
2688+ tmpVector);
2689+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, tmpVector);
2690+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
2691+ twoElec,
2692+ twiceMoJK,
2693+ tmpVector);
2694+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, tmpVector);
2695+ MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
2696+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
2697+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
2698+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
2699+ MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
2700+ MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
2701+ MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
2702+ MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
2703+ // End of algorithm using blas
2704+ */
2705+
2706+ /*
2707+ // Second algorithm using blas.
2708+ // This algorithm uses DGEMM.
2709+ double** twoElec = NULL;
2710+ double* twiceMoIJ = NULL;
2711+ double* twiceMoIK = NULL;
2712+ double* twiceMoIL = NULL;
2713+ double** twiceMoB = NULL;
2714+ double** tmpMatrix = NULL;
2715+ int numAOs = this->molecule->GetTotalNumberAOs();
2716+ MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
2717+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
2718+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
2719+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
2720+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
2721+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
2722+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
2723+ const Atom& atomA = *this->molecule->GetAtom(A);
2724+ int firstAOIndexA = atomA.GetFirstAOIndex();
2725+ int lastAOIndexA = atomA.GetLastAOIndex();
2726+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2727+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2728+ twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
2729+ twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
2730+ twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
2731+ }
2732+ }
2733+ }
2734+
2735+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
2736+ const Atom& atomB = *this->molecule->GetAtom(B);
2737+ int firstAOIndexB = atomB.GetFirstAOIndex();
2738+ int lastAOIndexB = atomB.GetLastAOIndex();
2739+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2740+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2741+ twiceMoB[0][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
2742+ twiceMoB[1][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
2743+ twiceMoB[2][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
2744+ }
2745+ }
2746+ }
2747+
2748+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
2749+ const Atom& atomA = *this->molecule->GetAtom(A);
2750+ int firstAOIndexA = atomA.GetFirstAOIndex();
2751+ int lastAOIndexA = atomA.GetLastAOIndex();
2752+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
2753+ const Atom& atomB = *this->molecule->GetAtom(B);
2754+ int firstAOIndexB = atomB.GetFirstAOIndex();
2755+ int lastAOIndexB = atomB.GetLastAOIndex();
2756+ double gamma = 0.0;
2757+ if(A!=B){
2758+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2759+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2760+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2761+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2762+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
2763+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
2764+ this->twoElecTwoCore[A]
2765+ [B]
2766+ [mu-firstAOIndexA]
2767+ [nu-firstAOIndexA]
2768+ [lambda-firstAOIndexB]
2769+ [sigma-firstAOIndexB];
2770+ }
2771+ }
2772+ }
2773+ }
2774+ }
2775+ else{
2776+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2777+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2778+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2779+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2780+ if(mu==nu && lambda==sigma){
2781+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2782+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
2783+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
2784+ }
2785+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
2786+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2787+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
2788+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
2789+ }
2790+ else{
2791+ gamma = 0.0;
2792+ }
2793+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
2794+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
2795+ }
2796+ }
2797+ }
2798+ }
2799+ }
2800+ }
2801+ }
2802+
2803+ MolDS_wrappers::Blas::GetInstance()->Dgemm(false, true, true,
2804+ this->molecule->GetNumberAtoms()*dxy*dxy,
2805+ 3,
2806+ this->molecule->GetNumberAtoms()*dxy*dxy,
2807+ 1.0,
2808+ twoElec,
2809+ twiceMoB,
2810+ 0.0,
2811+ tmpMatrix);
2812+ value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
2813+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
2814+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
2815+ MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
2816+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
2817+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
2818+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
2819+ MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
2820+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
2821+ // End of second algorithm using blas
2822+ */
2823+
2824+ /*
2825+ // slow algorithm
2826+ value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL,
2827+ *this->molecule,
2828+ this->fockMatrix, NULL)
2829+ -1.0*this->GetMolecularIntegralElement(moI, moK, moJ, moL,
2830+ *this->molecule,
2831+ this->fockMatrix, NULL)
2832+ -1.0*this->GetMolecularIntegralElement(moI, moL, moJ, moK,
2833+ *this->molecule,
2834+ this->fockMatrix, NULL);
2835+ */
2836+ return value;
2837+}
2838+
24372839 void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
24382840 const Molecule& molecule) const{
24392841 #ifdef MOLDS_DBG
@@ -6111,4 +6513,3 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
61116513
61126514 }
61136515
6114-
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -101,6 +101,7 @@ protected:
101101 double const* const* fockMatrix,
102102 double const* const* gammaAB) const;
103103 virtual void CalcCISMatrix(double** matrixCIS) const;
104+ virtual double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
104105 private:
105106 std::string errorMessageMultipoleA;
106107 std::string errorMessageMultipoleB;
--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -3417,9 +3417,12 @@ double ZindoS::GetKRElement(int moI, int moJ, int moK, int moL) const{
34173417 return 0.5*value;
34183418 }
34193419
3420+// see common term in eqs. (45) and (46) in [PT_1996],
3421+// that is, 4.0(ij|kl) - (ik|jl) - (il|jk).
34203422 double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
34213423 double value = 0.0;
34223424
3425+/*
34233426 // Fast algorith, but this is not easy to read.
34243427 // Slow algorithm is alos written below.
34253428 for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
@@ -3558,7 +3561,7 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
35583561 }
35593562 }
35603563 // End of the fast algorith.
3561-
3564+*/
35623565 /*
35633566 // Algorithm using blas
35643567 double** twoElec = NULL;
@@ -3802,7 +3805,6 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
38023805 // End of second algorithm using blas
38033806 */
38043807
3805- /*
38063808 // slow algorithm
38073809 value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL,
38083810 *this->molecule,
@@ -3813,7 +3815,6 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
38133815 -1.0*this->GetMolecularIntegralElement(moI, moL, moJ, moK,
38143816 *this->molecule,
38153817 this->fockMatrix, NULL);
3816- */
38173818 return value;
38183819 }
38193820
--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -116,7 +116,7 @@ protected:
116116 double GetKNRElement(int moI, int moJ, int moK, int moL) const;
117117 double GetKRElement(int moI, int moJ, int moK, int moL) const;
118118 double GetKRDagerElement(int moI, int moJ, int moK, int moL) const;
119- double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
119+ virtual double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
120120 void CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR,
121121 const std::vector<MoIndexPair>& nonRedundantQIndeces) const;
122122 void CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,