OSDN -- 開發和下載開源軟體
繁體中文翻譯狀態 translation status for zh-tw
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

提交

標籤
無標籤

Frequently used words (click to add to your profile)

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

Commit MetaInfo

修訂e3d5c38df4d7abc9f8d237c29d0a486dbeba90aa (tree)
時間2011-10-13 07:30:17
作者Mikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Mndo::GetNddoRepulsionIntegral is implemented. #26393

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

Change Summary

差異

--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -17,6 +17,7 @@ public:
1717 protected:
1818 virtual void SetMessages();
1919 virtual void SetEnableAtomTypes();
20+ virtual void CalcCoreRepulsionEnergy();
2021 virtual double GetFockDiagElement(Atom* atomA,
2122 int atomAIndex,
2223 int mu,
@@ -63,8 +64,9 @@ private:
6364 string errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles;
6465 string errorMessageMultipoleA;
6566 string errorMessageMultipoleB;
66- double GetNddoRepulsionIntegral(Atom* atomA, OrbitalType Mu, OrbitalType Nu,
67- Atom* atomB, OrbitalType Lambda, OrbitalType Sigma);
67+ string errorMessageGetNddoRepulsionIntegral;
68+ double GetNddoRepulsionIntegral(Atom* atomA, OrbitalType mu, OrbitalType nu,
69+ Atom* atomB, OrbitalType lambda, OrbitalType sigma);
6870 double GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
6971 MultipoleType multipoleB,
7072 double rhoA,
@@ -99,10 +101,12 @@ void Mndo::SetMessages(){
99101 this->errorMessageCalcCISMatrix
100102 = "Error in mndo::Mndo::CalcCISMatrix: Non available orbital is contained.\n";
101103 this->errorMessageDavidsonNotConverged = "Error in mndo::Mndo::DoesCISDavidson: Davidson did not met convergence criterion. \n";
102- this->errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles
104+ this->errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles
103105 = "Error in mndo:: Mndo::GetSemiEmpiricalMultipoleInteraction: Bad multipole combintaion is set\n";
104- this->errorMessageMultipoleA = "Multipole A is: ";
105- this->errorMessageMultipoleB = "Multipole B is: ";
106+ this->errorMessageMultipoleA = "Multipole A is: ";
107+ this->errorMessageMultipoleB = "Multipole B is: ";
108+ this->errorMessageGetNddoRepulsionIntegral = "Error in mndo::MNDO::GetNddoRepulsionIntegral: Bad orbital is set.\n";
109+
106110 this->messageSCFMetConvergence = "\n\n\n\t\tMNDO/S-SCF met convergence criterion(^^b\n\n\n";
107111 this->messageStartSCF = "********** START: MNDO/S-SCF **********\n";
108112 this->messageDoneSCF = "********** DONE: MNDO/S-SCF **********\n\n\n";
@@ -120,6 +124,41 @@ void Mndo::SetEnableAtomTypes(){
120124 this->enableAtomTypes.push_back(S);
121125 }
122126
127+void Mndo::CalcCoreRepulsionEnergy(){
128+ double energy = 0.0;
129+ double distance = 0.0;
130+ double twoElecInt = 0.0;
131+ double alphaA = 0.0;
132+ double alphaB = 0.0;
133+ Atom* atomA = NULL;
134+ Atom* atomB = NULL;
135+ double temp = 0.0;
136+ double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();
137+ for(int i=0; i<this->molecule->GetAtomVect()->size(); i++){
138+ atomA = (*this->molecule->GetAtomVect())[i];
139+ alphaA = atomA->GetMndoAlpha();
140+ for(int j=i+1; j<this->molecule->GetAtomVect()->size(); j++){
141+ atomB = (*this->molecule->GetAtomVect())[j];
142+ alphaB = atomB->GetMndoAlpha();
143+ distance = this->molecule->GetDistanceAtoms(i, j);
144+ if(atomA->GetAtomType() == H && (atomB->GetAtomType() == N ||
145+ atomB->GetAtomType() == O) ){
146+ temp = 1.0 + (distance/ang2AU)*exp(-alphaB*distance) + exp(-alphaA*distance);
147+ }
148+ else if(atomB->GetAtomType() == H && (atomA->GetAtomType() == N ||
149+ atomA->GetAtomType() == O) ){
150+ temp = 1.0 + (distance/ang2AU)*exp(-alphaA*distance) + exp(-alphaB*distance);
151+ }
152+ else{
153+ temp = 1.0 + exp(-alphaA*distance) + exp(-alphaB*distance);
154+ }
155+ twoElecInt = this->GetNddoRepulsionIntegral(atomA, s, s, atomB, s, s);
156+ energy += atomA->GetCoreCharge()*atomB->GetCoreCharge()*twoElecInt*temp;
157+ }
158+ }
159+ this->coreRepulsionEnergy = energy;
160+}
161+
123162 double Mndo::GetFockDiagElement(Atom* atomA, int atomAIndex, int mu,
124163 Molecule* molecule, double** gammaAB,
125164 double** orbitalElectronPopulation, double* atomicElectronPopulation,
@@ -481,11 +520,48 @@ double Mndo::GetElectronCoreAttraction(Atom* atomA, Atom* atomB,
481520 return -1.0*atomB->GetCoreCharge()*twoElecTwoCoreMatrixTwoAtoms[mu][nu][s][s];
482521 }
483522
523+void Mndo::CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB){
524+
525+ MolDS_cndo::Cndo2::CalcDiatomicOverlapInDiatomicFrame(diatomicOverlap, atomA, atomB);
526+
527+ /*
528+ for(int i=0;i<OrbitalType_end;i++){
529+ for(int j=0;j<OrbitalType_end;j++){
530+ printf("diatomicOverlap[%d][%d]=%lf\n",i,j,diatomicOverlap[i][j]);
531+ }
532+ }
533+ */
534+
535+}
536+
537+void Mndo::CalcDiatomicOverlapFirstDerivativeInDiatomicFrame(
538+ double** diatomicOverlapDeri,
539+ Atom* atomA, Atom* atomB){
540+
541+ MolDS_cndo::Cndo2::CalcDiatomicOverlapFirstDerivativeInDiatomicFrame(
542+ diatomicOverlapDeri,atomA, atomB);
543+
544+}
545+// The order of mol, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
546+double Mndo::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
547+ Molecule* molecule, double** fockMatrix, double** gammaAB){
548+ double value = 0.0;
549+ return value;
550+}
551+
552+void Mndo::CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir){
553+}
554+
555+// electronicStateIndex is index of the electroinc eigen state.
556+// "electronicStateIndex = 0" means electronic ground state.
557+void Mndo::CalcForce(int electronicStateIndex){
558+}
559+
484560 // See Apendix in [DT_1977]
485-// Orbital Mu and Nu belong atom A,
486-// orbital Lambda and Sigma belong atomB.
487-double Mndo::GetNddoRepulsionIntegral(Atom* atomA, OrbitalType Mu, OrbitalType Nu,
488- Atom* atomB, OrbitalType Lambda, OrbitalType Sigma){
561+// Orbital mu and nu belong atom A,
562+// orbital lambda and sigma belong atomB.
563+double Mndo::GetNddoRepulsionIntegral(Atom* atomA, OrbitalType mu, OrbitalType nu,
564+ Atom* atomB, OrbitalType lambda, OrbitalType sigma){
489565 double value = 0.0;
490566 double DA=0.0;
491567 double DB=0.0;
@@ -495,218 +571,619 @@ double Mndo::GetNddoRepulsionIntegral(Atom* atomA, OrbitalType Mu, OrbitalType N
495571 int lA = 0;
496572 int lB = 0;
497573 // (28) in [DT_1977]
498- if(Mu == s && Nu == s && Lambda == s && Sigma == s){
574+ if(mu == s && nu == s && lambda == s && sigma == s){
575+ DA = atomA->GetMndoDerivedParameterD(0);
576+ DB = atomB->GetMndoDerivedParameterD(0);
577+ rhoA = atomA->GetMndoDerivedParameterRho(0);
578+ rhoB = atomB->GetMndoDerivedParameterRho(0);
579+ value = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
499580 }
500581 // (29) in [DT_1977]
501- else if(Mu == s && Nu == s && Lambda == px && Sigma == px){
502- }
503- else if(Mu == s && Nu == s && Lambda == py && Sigma == py){
582+ else if(mu == s && nu == s && lambda == px && sigma == px){
583+ DA = atomA->GetMndoDerivedParameterD(0);
584+ DB = atomB->GetMndoDerivedParameterD(0);
585+ rhoA = atomA->GetMndoDerivedParameterRho(0);
586+ rhoB = atomB->GetMndoDerivedParameterRho(0);
587+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
588+ DA = atomA->GetMndoDerivedParameterD(0);
589+ DB = atomB->GetMndoDerivedParameterD(2);
590+ rhoA = atomA->GetMndoDerivedParameterRho(0);
591+ rhoB = atomB->GetMndoDerivedParameterRho(2);
592+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
593+ value = temp1 + temp2;
594+ }
595+ else if(mu == s && nu == s && lambda == py && sigma == py){
596+ DA = atomA->GetMndoDerivedParameterD(0);
597+ DB = atomB->GetMndoDerivedParameterD(0);
598+ rhoA = atomA->GetMndoDerivedParameterRho(0);
599+ rhoB = atomB->GetMndoDerivedParameterRho(0);
600+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
601+ DA = atomA->GetMndoDerivedParameterD(0);
602+ DB = atomB->GetMndoDerivedParameterD(2);
603+ rhoA = atomA->GetMndoDerivedParameterRho(0);
604+ rhoB = atomB->GetMndoDerivedParameterRho(2);
605+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
606+ value = temp1 + temp2;
504607 }
505608 // (30) in [DT_1977]
506- else if(Mu == s && Nu == s && Lambda == pz && Sigma == pz){
609+ else if(mu == s && nu == s && lambda == pz && sigma == pz){
610+ DA = atomA->GetMndoDerivedParameterD(0);
611+ DB = atomB->GetMndoDerivedParameterD(0);
612+ rhoA = atomA->GetMndoDerivedParameterRho(0);
613+ rhoB = atomB->GetMndoDerivedParameterRho(0);
614+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
615+ DA = atomA->GetMndoDerivedParameterD(0);
616+ DB = atomB->GetMndoDerivedParameterD(2);
617+ rhoA = atomA->GetMndoDerivedParameterRho(0);
618+ rhoB = atomB->GetMndoDerivedParameterRho(2);
619+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
620+ value = temp1 + temp2;
507621 }
508622 // (31) in [DT_1977]
509- else if(Mu == px && Nu == px && Lambda == s && Sigma == s){
510- }
511- else if(Mu == py && Nu == py && Lambda == s && Sigma == s){
623+ else if(mu == px && nu == px && lambda == s && sigma == s){
624+ DA = atomA->GetMndoDerivedParameterD(0);
625+ DB = atomB->GetMndoDerivedParameterD(0);
626+ rhoA = atomA->GetMndoDerivedParameterRho(0);
627+ rhoB = atomB->GetMndoDerivedParameterRho(0);
628+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
629+ DA = atomA->GetMndoDerivedParameterD(2);
630+ DB = atomB->GetMndoDerivedParameterD(0);
631+ rhoA = atomA->GetMndoDerivedParameterRho(2);
632+ rhoB = atomB->GetMndoDerivedParameterRho(0);
633+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
634+ value = temp1 + temp2;
635+ }
636+ else if(mu == py && nu == py && lambda == s && sigma == s){
637+ DA = atomA->GetMndoDerivedParameterD(0);
638+ DB = atomB->GetMndoDerivedParameterD(0);
639+ rhoA = atomA->GetMndoDerivedParameterRho(0);
640+ rhoB = atomB->GetMndoDerivedParameterRho(0);
641+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
642+ DA = atomA->GetMndoDerivedParameterD(2);
643+ DB = atomB->GetMndoDerivedParameterD(0);
644+ rhoA = atomA->GetMndoDerivedParameterRho(2);
645+ rhoB = atomB->GetMndoDerivedParameterRho(0);
646+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
647+ value = temp1 + temp2;
512648 }
513649 // (32) in [DT_1977]
514- else if(Mu == pz && Nu == pz && Lambda == s && Sigma == s){
650+ else if(mu == pz && nu == pz && lambda == s && sigma == s){
651+ DA = atomA->GetMndoDerivedParameterD(0);
652+ DB = atomB->GetMndoDerivedParameterD(0);
653+ rhoA = atomA->GetMndoDerivedParameterRho(0);
654+ rhoB = atomB->GetMndoDerivedParameterRho(0);
655+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
656+ DA = atomA->GetMndoDerivedParameterD(2);
657+ DB = atomB->GetMndoDerivedParameterD(0);
658+ rhoA = atomA->GetMndoDerivedParameterRho(2);
659+ rhoB = atomB->GetMndoDerivedParameterRho(0);
660+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
661+ value = temp1 + temp2;
515662 }
516663 // (33) in [DT_1977]
517- else if(Mu == px && Nu == px && Lambda == px && Sigma == px){
518- }
519- else if(Mu == py && Nu == py && Lambda == py && Sigma == py){
664+ else if(mu == px && nu == px && lambda == px && sigma == px){
665+ DA = atomA->GetMndoDerivedParameterD(0);
666+ DB = atomB->GetMndoDerivedParameterD(0);
667+ rhoA = atomA->GetMndoDerivedParameterRho(0);
668+ rhoB = atomB->GetMndoDerivedParameterRho(0);
669+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
670+ DA = atomA->GetMndoDerivedParameterD(0);
671+ DB = atomB->GetMndoDerivedParameterD(2);
672+ rhoA = atomA->GetMndoDerivedParameterRho(0);
673+ rhoB = atomB->GetMndoDerivedParameterRho(2);
674+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
675+ DA = atomA->GetMndoDerivedParameterD(2);
676+ DB = atomB->GetMndoDerivedParameterD(0);
677+ rhoA = atomA->GetMndoDerivedParameterRho(2);
678+ rhoB = atomB->GetMndoDerivedParameterRho(0);
679+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
680+ DA = atomA->GetMndoDerivedParameterD(2);
681+ DB = atomB->GetMndoDerivedParameterD(2);
682+ rhoA = atomA->GetMndoDerivedParameterRho(2);
683+ rhoB = atomB->GetMndoDerivedParameterRho(2);
684+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, Qxx, rhoA, rhoB, DA, DB, Rab);
685+ value = temp1 + temp2 + temp3 + temp4;
686+ }
687+ else if(mu == py && nu == py && lambda == py && sigma == py){
688+ DA = atomA->GetMndoDerivedParameterD(0);
689+ DB = atomB->GetMndoDerivedParameterD(0);
690+ rhoA = atomA->GetMndoDerivedParameterRho(0);
691+ rhoB = atomB->GetMndoDerivedParameterRho(0);
692+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
693+ DA = atomA->GetMndoDerivedParameterD(0);
694+ DB = atomB->GetMndoDerivedParameterD(2);
695+ rhoA = atomA->GetMndoDerivedParameterRho(0);
696+ rhoB = atomB->GetMndoDerivedParameterRho(2);
697+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
698+ DA = atomA->GetMndoDerivedParameterD(2);
699+ DB = atomB->GetMndoDerivedParameterD(0);
700+ rhoA = atomA->GetMndoDerivedParameterRho(2);
701+ rhoB = atomB->GetMndoDerivedParameterRho(0);
702+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
703+ DA = atomA->GetMndoDerivedParameterD(2);
704+ DB = atomB->GetMndoDerivedParameterD(2);
705+ rhoA = atomA->GetMndoDerivedParameterRho(2);
706+ rhoB = atomB->GetMndoDerivedParameterRho(2);
707+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, Qyy, rhoA, rhoB, DA, DB, Rab);
708+ value = temp1 + temp2 + temp3 + temp4;
520709 }
521710 // (34) in [DT_1977]
522- else if(Mu == px && Nu == px && Lambda == py && Sigma == py){
523- }
524- else if(Mu == py && Nu == py && Lambda == px && Sigma == px){
711+ else if(mu == px && nu == px && lambda == py && sigma == py){
712+ DA = atomA->GetMndoDerivedParameterD(0);
713+ DB = atomB->GetMndoDerivedParameterD(0);
714+ rhoA = atomA->GetMndoDerivedParameterRho(0);
715+ rhoB = atomB->GetMndoDerivedParameterRho(0);
716+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
717+ DA = atomA->GetMndoDerivedParameterD(0);
718+ DB = atomB->GetMndoDerivedParameterD(2);
719+ rhoA = atomA->GetMndoDerivedParameterRho(0);
720+ rhoB = atomB->GetMndoDerivedParameterRho(2);
721+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
722+ DA = atomA->GetMndoDerivedParameterD(2);
723+ DB = atomB->GetMndoDerivedParameterD(0);
724+ rhoA = atomA->GetMndoDerivedParameterRho(2);
725+ rhoB = atomB->GetMndoDerivedParameterRho(0);
726+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
727+ DA = atomA->GetMndoDerivedParameterD(2);
728+ DB = atomB->GetMndoDerivedParameterD(2);
729+ rhoA = atomA->GetMndoDerivedParameterRho(2);
730+ rhoB = atomB->GetMndoDerivedParameterRho(2);
731+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, Qyy, rhoA, rhoB, DA, DB, Rab);
732+ value = temp1 + temp2 + temp3 + temp4;
733+ }
734+ else if(mu == py && nu == py && lambda == px && sigma == px){
735+ DA = atomA->GetMndoDerivedParameterD(0);
736+ DB = atomB->GetMndoDerivedParameterD(0);
737+ rhoA = atomA->GetMndoDerivedParameterRho(0);
738+ rhoB = atomB->GetMndoDerivedParameterRho(0);
739+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
740+ DA = atomA->GetMndoDerivedParameterD(0);
741+ DB = atomB->GetMndoDerivedParameterD(2);
742+ rhoA = atomA->GetMndoDerivedParameterRho(0);
743+ rhoB = atomB->GetMndoDerivedParameterRho(2);
744+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
745+ DA = atomA->GetMndoDerivedParameterD(2);
746+ DB = atomB->GetMndoDerivedParameterD(0);
747+ rhoA = atomA->GetMndoDerivedParameterRho(2);
748+ rhoB = atomB->GetMndoDerivedParameterRho(0);
749+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
750+ DA = atomA->GetMndoDerivedParameterD(2);
751+ DB = atomB->GetMndoDerivedParameterD(2);
752+ rhoA = atomA->GetMndoDerivedParameterRho(2);
753+ rhoB = atomB->GetMndoDerivedParameterRho(2);
754+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, Qxx, rhoA, rhoB, DA, DB, Rab);
755+ value = temp1 + temp2 + temp3 + temp4;
525756 }
526757 // (35) in [DT_1977]
527- else if(Mu == px && Nu == px && Lambda == pz && Sigma == pz){
528- }
529- else if(Mu == py && Nu == py && Lambda == pz && Sigma == pz){
758+ else if(mu == px && nu == px && lambda == pz && sigma == pz){
759+ DA = atomA->GetMndoDerivedParameterD(0);
760+ DB = atomB->GetMndoDerivedParameterD(0);
761+ rhoA = atomA->GetMndoDerivedParameterRho(0);
762+ rhoB = atomB->GetMndoDerivedParameterRho(0);
763+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
764+ DA = atomA->GetMndoDerivedParameterD(0);
765+ DB = atomB->GetMndoDerivedParameterD(2);
766+ rhoA = atomA->GetMndoDerivedParameterRho(0);
767+ rhoB = atomB->GetMndoDerivedParameterRho(2);
768+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
769+ DA = atomA->GetMndoDerivedParameterD(2);
770+ DB = atomB->GetMndoDerivedParameterD(0);
771+ rhoA = atomA->GetMndoDerivedParameterRho(2);
772+ rhoB = atomB->GetMndoDerivedParameterRho(0);
773+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
774+ DA = atomA->GetMndoDerivedParameterD(2);
775+ DB = atomB->GetMndoDerivedParameterD(2);
776+ rhoA = atomA->GetMndoDerivedParameterRho(2);
777+ rhoB = atomB->GetMndoDerivedParameterRho(2);
778+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, Qzz, rhoA, rhoB, DA, DB, Rab);
779+ value = temp1 + temp2 + temp3 + temp4;
780+ }
781+ else if(mu == py && nu == py && lambda == pz && sigma == pz){
782+ DA = atomA->GetMndoDerivedParameterD(0);
783+ DB = atomB->GetMndoDerivedParameterD(0);
784+ rhoA = atomA->GetMndoDerivedParameterRho(0);
785+ rhoB = atomB->GetMndoDerivedParameterRho(0);
786+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
787+ DA = atomA->GetMndoDerivedParameterD(0);
788+ DB = atomB->GetMndoDerivedParameterD(2);
789+ rhoA = atomA->GetMndoDerivedParameterRho(0);
790+ rhoB = atomB->GetMndoDerivedParameterRho(2);
791+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
792+ DA = atomA->GetMndoDerivedParameterD(2);
793+ DB = atomB->GetMndoDerivedParameterD(0);
794+ rhoA = atomA->GetMndoDerivedParameterRho(2);
795+ rhoB = atomB->GetMndoDerivedParameterRho(0);
796+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
797+ DA = atomA->GetMndoDerivedParameterD(2);
798+ DB = atomB->GetMndoDerivedParameterD(2);
799+ rhoA = atomA->GetMndoDerivedParameterRho(2);
800+ rhoB = atomB->GetMndoDerivedParameterRho(2);
801+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, Qzz, rhoA, rhoB, DA, DB, Rab);
802+ value = temp1 + temp2 + temp3 + temp4;
530803 }
531804 // (36) in [DT_1977]
532- else if(Mu == pz && Nu == pz && Lambda == px && Sigma == px){
533- }
534- else if(Mu == pz && Nu == pz && Lambda == py && Sigma == py){
805+ else if(mu == pz && nu == pz && lambda == px && sigma == px){
806+ DA = atomA->GetMndoDerivedParameterD(0);
807+ DB = atomB->GetMndoDerivedParameterD(0);
808+ rhoA = atomA->GetMndoDerivedParameterRho(0);
809+ rhoB = atomB->GetMndoDerivedParameterRho(0);
810+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
811+ DA = atomA->GetMndoDerivedParameterD(0);
812+ DB = atomB->GetMndoDerivedParameterD(2);
813+ rhoA = atomA->GetMndoDerivedParameterRho(0);
814+ rhoB = atomB->GetMndoDerivedParameterRho(2);
815+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
816+ DA = atomA->GetMndoDerivedParameterD(2);
817+ DB = atomB->GetMndoDerivedParameterD(0);
818+ rhoA = atomA->GetMndoDerivedParameterRho(2);
819+ rhoB = atomB->GetMndoDerivedParameterRho(0);
820+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
821+ DA = atomA->GetMndoDerivedParameterD(2);
822+ DB = atomB->GetMndoDerivedParameterD(2);
823+ rhoA = atomA->GetMndoDerivedParameterRho(2);
824+ rhoB = atomB->GetMndoDerivedParameterRho(2);
825+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, Qxx, rhoA, rhoB, DA, DB, Rab);
826+ value = temp1 + temp2 + temp3 + temp4;
827+ }
828+ else if(mu == pz && nu == pz && lambda == py && sigma == py){
829+ DA = atomA->GetMndoDerivedParameterD(0);
830+ DB = atomB->GetMndoDerivedParameterD(0);
831+ rhoA = atomA->GetMndoDerivedParameterRho(0);
832+ rhoB = atomB->GetMndoDerivedParameterRho(0);
833+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
834+ DA = atomA->GetMndoDerivedParameterD(0);
835+ DB = atomB->GetMndoDerivedParameterD(2);
836+ rhoA = atomA->GetMndoDerivedParameterRho(0);
837+ rhoB = atomB->GetMndoDerivedParameterRho(2);
838+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
839+ DA = atomA->GetMndoDerivedParameterD(2);
840+ DB = atomB->GetMndoDerivedParameterD(0);
841+ rhoA = atomA->GetMndoDerivedParameterRho(2);
842+ rhoB = atomB->GetMndoDerivedParameterRho(0);
843+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
844+ DA = atomA->GetMndoDerivedParameterD(2);
845+ DB = atomB->GetMndoDerivedParameterD(2);
846+ rhoA = atomA->GetMndoDerivedParameterRho(2);
847+ rhoB = atomB->GetMndoDerivedParameterRho(2);
848+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, Qyy, rhoA, rhoB, DA, DB, Rab);
849+ value = temp1 + temp2 + temp3 + temp4;
535850 }
536851 // (37) in [DT_1977]
537- else if(Mu == pz && Nu == pz && Lambda == pz && Sigma == pz){
852+ else if(mu == pz && nu == pz && lambda == pz && sigma == pz){
853+ DA = atomA->GetMndoDerivedParameterD(0);
854+ DB = atomB->GetMndoDerivedParameterD(0);
855+ rhoA = atomA->GetMndoDerivedParameterRho(0);
856+ rhoB = atomB->GetMndoDerivedParameterRho(0);
857+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, sQ, rhoA, rhoB, DA, DB, Rab);
858+ DA = atomA->GetMndoDerivedParameterD(0);
859+ DB = atomB->GetMndoDerivedParameterD(2);
860+ rhoA = atomA->GetMndoDerivedParameterRho(0);
861+ rhoB = atomB->GetMndoDerivedParameterRho(2);
862+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
863+ DA = atomA->GetMndoDerivedParameterD(2);
864+ DB = atomB->GetMndoDerivedParameterD(0);
865+ rhoA = atomA->GetMndoDerivedParameterRho(2);
866+ rhoB = atomB->GetMndoDerivedParameterRho(0);
867+ double temp3 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
868+ DA = atomA->GetMndoDerivedParameterD(2);
869+ DB = atomB->GetMndoDerivedParameterD(2);
870+ rhoA = atomA->GetMndoDerivedParameterRho(2);
871+ rhoB = atomB->GetMndoDerivedParameterRho(2);
872+ double temp4 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, Qzz, rhoA, rhoB, DA, DB, Rab);
873+ value = temp1 + temp2 + temp3 + temp4;
538874 }
539875 // (38) in [DT_1977]
540- else if(Mu == s && Nu == pz && Lambda == s && Sigma == s){
876+ else if(mu == s && nu == pz && lambda == s && sigma == s){
877+ DA = atomA->GetMndoDerivedParameterD(1);
878+ DB = atomB->GetMndoDerivedParameterD(0);
879+ rhoA = atomA->GetMndoDerivedParameterRho(1);
880+ rhoB = atomB->GetMndoDerivedParameterRho(0);
881+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muz, sQ, rhoA, rhoB, DA, DB, Rab);
882+ value = temp1;
541883 }
542- else if(Mu == pz && Nu == s && Lambda == s && Sigma == s){
543- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
884+ else if(mu == pz && nu == s && lambda == s && sigma == s){
885+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
544886 }
545887 // (39) in [DT_1977]
546- else if(Mu == s && Nu == pz && Lambda == px && Sigma == px){
547- }
548- else if(Mu == pz && Nu == s && Lambda == px && Sigma == px){
549- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
550- }
551- else if(Mu == s && Nu == pz && Lambda == py && Sigma == py){
552- }
553- else if(Mu == pz && Nu == s && Lambda == py && Sigma == py){
554- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
888+ else if(mu == s && nu == pz && lambda == px && sigma == px){
889+ DA = atomA->GetMndoDerivedParameterD(1);
890+ DB = atomB->GetMndoDerivedParameterD(0);
891+ rhoA = atomA->GetMndoDerivedParameterRho(1);
892+ rhoB = atomB->GetMndoDerivedParameterRho(0);
893+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muz, sQ, rhoA, rhoB, DA, DB, Rab);
894+ DA = atomA->GetMndoDerivedParameterD(1);
895+ DB = atomB->GetMndoDerivedParameterD(2);
896+ rhoA = atomA->GetMndoDerivedParameterRho(1);
897+ rhoB = atomB->GetMndoDerivedParameterRho(2);
898+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(muz, Qxx, rhoA, rhoB, DA, DB, Rab);
899+ value = temp1 + temp2;
900+ }
901+ else if(mu == pz && nu == s && lambda == px && sigma == px){
902+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
903+ }
904+ else if(mu == s && nu == pz && lambda == py && sigma == py){
905+ DA = atomA->GetMndoDerivedParameterD(1);
906+ DB = atomB->GetMndoDerivedParameterD(0);
907+ rhoA = atomA->GetMndoDerivedParameterRho(1);
908+ rhoB = atomB->GetMndoDerivedParameterRho(0);
909+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muz, sQ, rhoA, rhoB, DA, DB, Rab);
910+ DA = atomA->GetMndoDerivedParameterD(1);
911+ DB = atomB->GetMndoDerivedParameterD(2);
912+ rhoA = atomA->GetMndoDerivedParameterRho(1);
913+ rhoB = atomB->GetMndoDerivedParameterRho(2);
914+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(muz, Qyy, rhoA, rhoB, DA, DB, Rab);
915+ value = temp1 + temp2;
916+ }
917+ else if(mu == pz && nu == s && lambda == py && sigma == py){
918+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
555919 }
556920 // (40) in [DT_1977]
557- else if(Mu == s && Nu == pz && Lambda == pz && Sigma == pz){
558- }
559- else if(Mu == pz && Nu == s && Lambda == pz && Sigma == pz){
560- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
921+ else if(mu == s && nu == pz && lambda == pz && sigma == pz){
922+ DA = atomA->GetMndoDerivedParameterD(1);
923+ DB = atomB->GetMndoDerivedParameterD(0);
924+ rhoA = atomA->GetMndoDerivedParameterRho(1);
925+ rhoB = atomB->GetMndoDerivedParameterRho(0);
926+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muz, sQ, rhoA, rhoB, DA, DB, Rab);
927+ DA = atomA->GetMndoDerivedParameterD(1);
928+ DB = atomB->GetMndoDerivedParameterD(2);
929+ rhoA = atomA->GetMndoDerivedParameterRho(1);
930+ rhoB = atomB->GetMndoDerivedParameterRho(2);
931+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(muz, Qzz, rhoA, rhoB, DA, DB, Rab);
932+ value = temp1 + temp2;
933+ }
934+ else if(mu == pz && nu == s && lambda == pz && sigma == pz){
935+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
561936 }
562937 // (41) in [DT_1977]
563- else if(Mu == s && Nu == s && Lambda == s && Sigma == pz){
938+ else if(mu == s && nu == s && lambda == s && sigma == pz){
939+ DA = atomA->GetMndoDerivedParameterD(0);
940+ DB = atomB->GetMndoDerivedParameterD(1);
941+ rhoA = atomA->GetMndoDerivedParameterRho(0);
942+ rhoB = atomB->GetMndoDerivedParameterRho(1);
943+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, muz, rhoA, rhoB, DA, DB, Rab);
944+ value = temp1;
564945 }
565- else if(Mu == s && Nu == s && Lambda == pz && Sigma == s){
566- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
946+ else if(mu == s && nu == s && lambda == pz && sigma == s){
947+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
567948 }
568949 // (42) in [DT_1977]
569- else if(Mu == px && Nu == px && Lambda == s && Sigma == pz){
570- }
571- else if(Mu == px && Nu == px && Lambda == pz && Sigma == s){
572- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
573- }
574- else if(Mu == py && Nu == py && Lambda == s && Sigma == pz){
575- }
576- else if(Mu == py && Nu == py && Lambda == pz && Sigma == s){
577- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
950+ else if(mu == px && nu == px && lambda == s && sigma == pz){
951+ DA = atomA->GetMndoDerivedParameterD(0);
952+ DB = atomB->GetMndoDerivedParameterD(1);
953+ rhoA = atomA->GetMndoDerivedParameterRho(0);
954+ rhoB = atomB->GetMndoDerivedParameterRho(1);
955+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, muz, rhoA, rhoB, DA, DB, Rab);
956+ DA = atomA->GetMndoDerivedParameterD(2);
957+ DB = atomB->GetMndoDerivedParameterD(1);
958+ rhoA = atomA->GetMndoDerivedParameterRho(2);
959+ rhoB = atomB->GetMndoDerivedParameterRho(1);
960+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qxx, muz, rhoA, rhoB, DA, DB, Rab);
961+ value = temp1 + temp2;
962+ }
963+ else if(mu == px && nu == px && lambda == pz && sigma == s){
964+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
965+ }
966+ else if(mu == py && nu == py && lambda == s && sigma == pz){
967+ DA = atomA->GetMndoDerivedParameterD(0);
968+ DB = atomB->GetMndoDerivedParameterD(1);
969+ rhoA = atomA->GetMndoDerivedParameterRho(0);
970+ rhoB = atomB->GetMndoDerivedParameterRho(1);
971+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, muz, rhoA, rhoB, DA, DB, Rab);
972+ DA = atomA->GetMndoDerivedParameterD(2);
973+ DB = atomB->GetMndoDerivedParameterD(1);
974+ rhoA = atomA->GetMndoDerivedParameterRho(2);
975+ rhoB = atomB->GetMndoDerivedParameterRho(1);
976+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qyy, muz, rhoA, rhoB, DA, DB, Rab);
977+ value = temp1 + temp2;
978+ }
979+ else if(mu == py && nu == py && lambda == pz && sigma == s){
980+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
578981 }
579982 // (43) in [DT_1977]
580- else if(Mu == pz && Nu == pz && Lambda == s && Sigma == pz){
581- }
582- else if(Mu == pz && Nu == pz && Lambda == pz && Sigma == s){
583- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
983+ else if(mu == pz && nu == pz && lambda == s && sigma == pz){
984+ DA = atomA->GetMndoDerivedParameterD(0);
985+ DB = atomB->GetMndoDerivedParameterD(1);
986+ rhoA = atomA->GetMndoDerivedParameterRho(0);
987+ rhoB = atomB->GetMndoDerivedParameterRho(1);
988+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(sQ, muz, rhoA, rhoB, DA, DB, Rab);
989+ DA = atomA->GetMndoDerivedParameterD(2);
990+ DB = atomB->GetMndoDerivedParameterD(1);
991+ rhoA = atomA->GetMndoDerivedParameterRho(2);
992+ rhoB = atomB->GetMndoDerivedParameterRho(1);
993+ double temp2 = this->GetSemiEmpiricalMultipoleInteraction(Qzz, muz, rhoA, rhoB, DA, DB, Rab);
994+ value = temp1 + temp2;
995+ }
996+ else if(mu == pz && nu == pz && lambda == pz && sigma == s){
997+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
584998 }
585999 // (44) in [DT_1977]
586- else if(Mu == s && Nu == px && Lambda == s && Sigma == px){
1000+ else if(mu == s && nu == px && lambda == s && sigma == px){
1001+ DA = atomA->GetMndoDerivedParameterD(1);
1002+ DB = atomB->GetMndoDerivedParameterD(1);
1003+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1004+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1005+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(mux, mux, rhoA, rhoB, DA, DB, Rab);
1006+ value = temp1;
5871007 }
588- else if(Mu == px && Nu == s && Lambda == s && Sigma == px){
589- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1008+ else if(mu == px && nu == s && lambda == s && sigma == px){
1009+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
5901010 }
591- else if(Mu == s && Nu == px && Lambda == px && Sigma == s){
592- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1011+ else if(mu == s && nu == px && lambda == px && sigma == s){
1012+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
5931013 }
594- else if(Mu == px && Nu == s && Lambda == px && Sigma == s){
595- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1014+ else if(mu == px && nu == s && lambda == px && sigma == s){
1015+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
5961016 }
597- else if(Mu == s && Nu == py && Lambda == s && Sigma == py){
1017+ else if(mu == s && nu == py && lambda == s && sigma == py){
1018+ DA = atomA->GetMndoDerivedParameterD(1);
1019+ DB = atomB->GetMndoDerivedParameterD(1);
1020+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1021+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1022+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muy, muy, rhoA, rhoB, DA, DB, Rab);
1023+ value = temp1;
5981024 }
599- else if(Mu == py && Nu == s && Lambda == s && Sigma == py){
600- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1025+ else if(mu == py && nu == s && lambda == s && sigma == py){
1026+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6011027 }
602- else if(Mu == s && Nu == py && Lambda == py && Sigma == s){
603- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1028+ else if(mu == s && nu == py && lambda == py && sigma == s){
1029+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6041030 }
605- else if(Mu == py && Nu == s && Lambda == py && Sigma == s){
606- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1031+ else if(mu == py && nu == s && lambda == py && sigma == s){
1032+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6071033 }
6081034 // (45) in [DT_1977]
609- else if(Mu == s && Nu == pz && Lambda == s && Sigma == pz){
1035+ else if(mu == s && nu == pz && lambda == s && sigma == pz){
1036+ DA = atomA->GetMndoDerivedParameterD(1);
1037+ DB = atomB->GetMndoDerivedParameterD(1);
1038+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1039+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1040+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muz, muz, rhoA, rhoB, DA, DB, Rab);
1041+ value = temp1;
6101042 }
611- else if(Mu == pz && Nu == s && Lambda == s && Sigma == pz){
612- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1043+ else if(mu == pz && nu == s && lambda == s && sigma == pz){
1044+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6131045 }
614- else if(Mu == s && Nu == pz && Lambda == pz && Sigma == s){
615- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1046+ else if(mu == s && nu == pz && lambda == pz && sigma == s){
1047+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6161048 }
617- else if(Mu == pz && Nu == s && Lambda == pz && Sigma == s){
618- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1049+ else if(mu == pz && nu == s && lambda == pz && sigma == s){
1050+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6191051 }
6201052 // (46) in [DT_1977]
621- else if(Mu == s && Nu == px && Lambda == px && Sigma == pz){
1053+ else if(mu == s && nu == px && lambda == px && sigma == pz){
1054+ DA = atomA->GetMndoDerivedParameterD(1);
1055+ DB = atomB->GetMndoDerivedParameterD(2);
1056+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1057+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1058+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(mux, Qxz, rhoA, rhoB, DA, DB, Rab);
1059+ value = temp1;
6221060 }
623- else if(Mu == px && Nu == s && Lambda == px && Sigma == pz){
624- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1061+ else if(mu == px && nu == s && lambda == px && sigma == pz){
1062+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6251063 }
626- else if(Mu == s && Nu == px && Lambda == pz && Sigma == px){
627- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1064+ else if(mu == s && nu == px && lambda == pz && sigma == px){
1065+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6281066 }
629- else if(Mu == px && Nu == s && Lambda == pz && Sigma == px){
630- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1067+ else if(mu == px && nu == s && lambda == pz && sigma == px){
1068+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6311069 }
632- else if(Mu == s && Nu == py && Lambda == py && Sigma == pz){
1070+ else if(mu == s && nu == py && lambda == py && sigma == pz){
1071+ DA = atomA->GetMndoDerivedParameterD(1);
1072+ DB = atomB->GetMndoDerivedParameterD(2);
1073+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1074+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1075+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(muy, Qyz, rhoA, rhoB, DA, DB, Rab);
1076+ value = temp1;
6331077 }
634- else if(Mu == py && Nu == s && Lambda == py && Sigma == pz){
635- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1078+ else if(mu == py && nu == s && lambda == py && sigma == pz){
1079+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6361080 }
637- else if(Mu == s && Nu == py && Lambda == pz && Sigma == py){
638- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1081+ else if(mu == s && nu == py && lambda == pz && sigma == py){
1082+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6391083 }
640- else if(Mu == py && Nu == s && Lambda == pz && Sigma == py){
641- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1084+ else if(mu == py && nu == s && lambda == pz && sigma == py){
1085+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6421086 }
6431087 // (47) in [DT_1977]
644- else if(Mu == px && Nu == pz && Lambda == s && Sigma == px){
1088+ else if(mu == px && nu == pz && lambda == s && sigma == px){
1089+ DA = atomA->GetMndoDerivedParameterD(2);
1090+ DB = atomB->GetMndoDerivedParameterD(1);
1091+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1092+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1093+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(Qxz, mux, rhoA, rhoB, DA, DB, Rab);
1094+ value = temp1;
6451095 }
646- else if(Mu == pz && Nu == px && Lambda == s && Sigma == px){
647- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1096+ else if(mu == pz && nu == px && lambda == s && sigma == px){
1097+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6481098 }
649- else if(Mu == px && Nu == pz && Lambda == px && Sigma == s){
650- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1099+ else if(mu == px && nu == pz && lambda == px && sigma == s){
1100+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6511101 }
652- else if(Mu == pz && Nu == px && Lambda == px && Sigma == s){
653- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1102+ else if(mu == pz && nu == px && lambda == px && sigma == s){
1103+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6541104 }
655- else if(Mu == py && Nu == pz && Lambda == s && Sigma == py){
1105+ else if(mu == py && nu == pz && lambda == s && sigma == py){
1106+ DA = atomA->GetMndoDerivedParameterD(2);
1107+ DB = atomB->GetMndoDerivedParameterD(1);
1108+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1109+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1110+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(Qyz, muy, rhoA, rhoB, DA, DB, Rab);
1111+ value = temp1;
6561112 }
657- else if(Mu == pz && Nu == py && Lambda == s && Sigma == py){
658- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1113+ else if(mu == pz && nu == py && lambda == s && sigma == py){
1114+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6591115 }
660- else if(Mu == py && Nu == pz && Lambda == py && Sigma == s){
661- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1116+ else if(mu == py && nu == pz && lambda == py && sigma == s){
1117+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6621118 }
663- else if(Mu == pz && Nu == py && Lambda == py && Sigma == s){
664- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1119+ else if(mu == pz && nu == py && lambda == py && sigma == s){
1120+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6651121 }
6661122 // (48) in [DT_1977]
667- else if(Mu == px && Nu == pz && Lambda == px && Sigma == pz){
1123+ else if(mu == px && nu == pz && lambda == px && sigma == pz){
1124+ DA = atomA->GetMndoDerivedParameterD(2);
1125+ DB = atomB->GetMndoDerivedParameterD(2);
1126+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1127+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1128+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(Qxz, Qxz, rhoA, rhoB, DA, DB, Rab);
1129+ value = temp1;
6681130 }
669- else if(Mu == pz && Nu == px && Lambda == px && Sigma == pz){
670- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1131+ else if(mu == pz && nu == px && lambda == px && sigma == pz){
1132+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6711133 }
672- else if(Mu == px && Nu == pz && Lambda == pz && Sigma == px){
673- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1134+ else if(mu == px && nu == pz && lambda == pz && sigma == px){
1135+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6741136 }
675- else if(Mu == pz && Nu == px && Lambda == pz && Sigma == px){
676- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1137+ else if(mu == pz && nu == px && lambda == pz && sigma == px){
1138+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6771139 }
678- else if(Mu == py && Nu == pz && Lambda == py && Sigma == pz){
1140+ else if(mu == py && nu == pz && lambda == py && sigma == pz){
1141+ DA = atomA->GetMndoDerivedParameterD(2);
1142+ DB = atomB->GetMndoDerivedParameterD(2);
1143+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1144+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1145+ double temp1 = this->GetSemiEmpiricalMultipoleInteraction(Qyz, Qyz, rhoA, rhoB, DA, DB, Rab);
1146+ value = temp1;
6791147 }
680- else if(Mu == pz && Nu == py && Lambda == py && Sigma == pz){
681- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1148+ else if(mu == pz && nu == py && lambda == py && sigma == pz){
1149+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6821150 }
683- else if(Mu == py && Nu == pz && Lambda == pz && Sigma == py){
684- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1151+ else if(mu == py && nu == pz && lambda == pz && sigma == py){
1152+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6851153 }
686- else if(Mu == pz && Nu == py && Lambda == pz && Sigma == py){
687- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1154+ else if(mu == pz && nu == py && lambda == pz && sigma == py){
1155+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
6881156 }
6891157 // (49) in [DT_1977] and p19 in [MOPAC_1990]
690- else if(Mu == px && Nu == py && Lambda == px && Sigma == py){
691- value = 0.5*(this->GetNddoRepulsionIntegral(atomA, Mu, Mu, atomB, Mu, Mu)
692- -this->GetNddoRepulsionIntegral(atomA, Mu, Mu, atomB, Nu, Nu));
1158+ else if(mu == px && nu == py && lambda == px && sigma == py){
1159+ value = 0.5*(this->GetNddoRepulsionIntegral(atomA, mu, mu, atomB, mu, mu)
1160+ -this->GetNddoRepulsionIntegral(atomA, mu, mu, atomB, nu, nu));
6931161 }
694- else if(Mu == py && Nu == px && Lambda == px && Sigma == py){
695- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Lambda, Sigma);
1162+ else if(mu == py && nu == px && lambda == px && sigma == py){
1163+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, lambda, sigma);
6961164 }
697- else if(Mu == px && Nu == py && Lambda == py && Sigma == px){
698- value = this->GetNddoRepulsionIntegral(atomA, Mu, Nu, atomB, Sigma, Lambda);
1165+ else if(mu == px && nu == py && lambda == py && sigma == px){
1166+ value = this->GetNddoRepulsionIntegral(atomA, mu, nu, atomB, sigma, lambda);
6991167 }
700- else if(Mu == py && Nu == px && Lambda == py && Sigma == px){
701- value = this->GetNddoRepulsionIntegral(atomA, Nu, Mu, atomB, Sigma, Lambda);
1168+ else if(mu == py && nu == px && lambda == py && sigma == px){
1169+ value = this->GetNddoRepulsionIntegral(atomA, nu, mu, atomB, sigma, lambda);
7021170 }
7031171 // d-orbitals
704- else if(Mu == dxy || Mu == dyz || Mu == dzz || Mu == dzx || Mu == dxxyy ||
705- Nu == dxy || Nu == dyz || Nu == dzz || Nu == dzx || Nu == dxxyy ||
706- Lambda == dxy || Lambda == dyz || Lambda == dzz || Lambda == dzx || Lambda == dxxyy ||
707- Sigma == dxy || Sigma == dyz || Sigma == dzz || Sigma == dzx || Sigma == dxxyy){
1172+ else if(mu == dxy || mu == dyz || mu == dzz || mu == dzx || mu == dxxyy ||
1173+ nu == dxy || nu == dyz || nu == dzz || nu == dzx || nu == dxxyy ||
1174+ lambda == dxy || lambda == dyz || lambda == dzz || lambda == dzx || lambda == dxxyy ||
1175+ sigma == dxy || sigma == dyz || sigma == dzz || sigma == dzx || sigma == dxxyy){
7081176
7091177 // ToDo: error log.
1178+ stringstream ss;
1179+ ss << this->errorMessageGetNddoRepulsionIntegral;
1180+ ss << this->errorMessageAtomA << AtomTypeStr(atomA->GetAtomType()) << endl;
1181+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(mu) << endl;
1182+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(nu) << endl;
1183+ ss << this->errorMessageAtomB << AtomTypeStr(atomB->GetAtomType()) << endl;
1184+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(lambda) << endl;
1185+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(sigma) << endl;
1186+ throw MolDSException(ss.str());
7101187 }
7111188 else{
7121189 value = 0.0;
@@ -948,43 +1425,6 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
9481425 return value;
9491426 }
9501427
951-void Mndo::CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB){
952-
953- MolDS_cndo::Cndo2::CalcDiatomicOverlapInDiatomicFrame(diatomicOverlap, atomA, atomB);
954-
955- /*
956- for(int i=0;i<OrbitalType_end;i++){
957- for(int j=0;j<OrbitalType_end;j++){
958- printf("diatomicOverlap[%d][%d]=%lf\n",i,j,diatomicOverlap[i][j]);
959- }
960- }
961- */
962-
963-}
964-
965-void Mndo::CalcDiatomicOverlapFirstDerivativeInDiatomicFrame(
966- double** diatomicOverlapDeri,
967- Atom* atomA, Atom* atomB){
968-
969- MolDS_cndo::Cndo2::CalcDiatomicOverlapFirstDerivativeInDiatomicFrame(
970- diatomicOverlapDeri,atomA, atomB);
971-
972-}
973-// The order of mol, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
974-double Mndo::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
975- Molecule* molecule, double** fockMatrix, double** gammaAB){
976- double value = 0.0;
977- return value;
978-}
979-
980-void Mndo::CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir){
981-}
982-
983-// electronicStateIndex is index of the electroinc eigen state.
984-// "electronicStateIndex = 0" means electronic ground state.
985-void Mndo::CalcForce(int electronicStateIndex){
986-}
987-
9881428 }
9891429 #endif
9901430