コミット: 7ce00af20ec2a294cf61e022258c3534bf666d50 - MolDS (git)

コミットメタ情報

• リビジョン: 7ce00af20ec2a294cf61e022258c3534bf666d50 (tree)
• 日時: 2010-12-20 19:12:56
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• コミッター: Mikiya Fujii

ログメッセージ

ZINDO/S is implemented, although it is not completed yet.

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

変更サマリ

• modified: doc/refferences.txt
• modified: src/base/atoms/Atom.h
• modified: src/base/atoms/Catom.h
• modified: src/base/atoms/Hatom.h
• modified: src/base/atoms/Liatom.h
• modified: src/base/atoms/Satom.h
• modified: src/cndo/Cndo2.h
• modified: src/indo/Indo.h
• modified: src/input.in
• modified: src/zindo/ZindoS.h

差分

--- a/doc/refferences.txt
+++ b/doc/refferences.txt
@@ -1,3 +1,4 @@
+[MN_1957] N. Mataga and K. Nishinoto, Z. Phys. Chem., 13, 140 (1957)
 [PSS_1965] J. A. Pople, D. P. Santry, and G. A. Segal, J. Chem. Phys., 43, S129 (1965)
 [PS_1965] J. A. Pople and G. A. Segal, J. Chem. Phys., 43, S136 (1965)
 [PS_1966] J. A. Pople and G. A. Segal, J. Chem. Phys., 44, 3289 (1966)
@@ -5,9 +6,9 @@
 [GD_1972] G. P.-Guillouzo and D. G. ET J. Deschamps, J. Mol. Strct., 14, 81 (1972)
 [RZ_1976] J. E. Ridley and M. C. Zerner, Theort. Chim. Acta (Berl.) 42, 223 (1976)
 [BZ_1979] A. D. Bacon and M. C. Zerner, Theort. Chim. Acta (Berl.) 53, 21 (1979)
+[HKLWNZ_1982] Z. S. Herman, R. F. Kirachner, G. H. Loew, U. T. M.-WesterHoff, A. Nazzal, and M. C. Zerner, Inorg. Chem., 21, 46 (1982)
 [AEZ_1986] W. P. Anderson, W. D. Edwards, and M. C. Zerner, Inorg. Chem. 25, 2728 (1986)
 [BFB_1997] M. A. Blanco, M. Fl{\'o}rez, and M. Bermejo, J. Mol. Struct. 419, 19 (1997)
 
 
 
-

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -41,6 +41,8 @@ protected:
    double imuAmuP;                      // Table 3.4 or 3.5 in J. A. Pople book
    double imuAmuD;                      // Table 3.4 or 3.5 in J. A. Pople book
    double bondingParameter;             // Table 3.2 and 3.4 in J. A. Pople book
+   double bondingParameterSZindo;        // Table 1 in [RZ_1976], table 1 in [HKLWNZ_1982], or table 3 in [AEZ_1986]
+   double bondingParameterDZindo;        // Table 1 in [RZ_1976], table 1 in [HKLWNZ_1982], or table 3 in [AEZ_1986]
    double coreCharge;                   // = Z_A in J. A. Pople book.
    double effectiveNuclearChargeK;      // Table 1.5 in J. A. Pople book
    double effectiveNuclearChargeL;      // Table 1.5 in J. A. Pople book
@@ -74,6 +76,7 @@ public:
    void SetXyz(double x, double y, double z);
    vector<OrbitalType> GetValence();
    double GetBondingParameter();
+   double GetBondingParameter(TheoryType theory, OrbitalType orbital);
    double GetCoreCharge();
    int GetFirstAOIndex();
    void SetFirstAOIndex(int firstAOIndex);
@@ -82,8 +85,31 @@ public:
    double GetImuAmu(OrbitalType orbitalType);  // return 0.5*(I_mu + A_mu)
    double GetOrbitalExponent(ShellType shellType, OrbitalType orbitalType);  // (1.73) in J. A. Pople book.
    virtual double GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory) = 0; // P82 - 83 in J. A. Pople book for INDO or Eq. (13) in [BZ_1979] for ZINDO/S
+   double GetCoreIntegral(OrbitalType orbital, bool isGuess, TheoryType theory);
    double GetIndoF2();
    double GetIndoG1();
+   double GetZindoF0ss();                // Table 1 in ref. [RZ_1976], Table 1 in [AEZ_1986], or Table 1 in [GD_1972]
+   double GetZindoF0sd();                  // Table 1 in [AEZ_1986]
+   double GetZindoF0dd();                  // Table 1 in [AEZ_1986]
+   double GetZindoG1sp();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoF2pp();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoG2sd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoG1pd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoF2pd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoG3pd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoF2dd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoF4dd();                 // Table 3 in ref. [BZ_1979]
+   double GetZindoF0ssLower();                 // Apendix in ref. [BZ_1979] 
+   double GetZindoF0sdLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoF0ddLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoG1spLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoF2ppLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoG2sdLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoG1pdLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoF2pdLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoG3pdLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoF2ddLower();                 // Apendix in ref. [BZ_1979]
+   double GetZindoF4ddLower();                 // Apendix in ref. [BZ_1979]
 };
 
 Atom::Atom(){
@@ -137,8 +163,32 @@ vector<OrbitalType> Atom::GetValence(){
    return this->valence;
 }
 
+double Atom::GetBondingParameter(TheoryType theory, OrbitalType orbital){
+
+   double value = 0.0;
+   if(theory == CNDO2 || theory == INDO){
+      value = this->bondingParameter;
+   }     
+   else if(theory == ZINDOS && ( orbital == s ||
+                                 orbital == px ||
+                                 orbital == py ||
+                                 orbital == pz ) ){
+      value = this->bondingParameterSZindo;
+   }
+   else if(theory == ZINDOS && ( orbital == dxy ||
+                                 orbital == dyz ||
+                                 orbital == dzz ||
+                                 orbital == dzx ||
+                                 orbital == dxxyy ) ){
+      value = this->bondingParameterDZindo;
+   }
+
+   return value;
+
+}
+
 double Atom::GetBondingParameter(){
-   return this->bondingParameter;
+   return this->GetBondingParameter(CNDO2, s);
 }
 
 double Atom::GetCoreCharge(){
@@ -268,7 +318,7 @@ double Atom::GetJpd(){
 
 // Part of Eq. (13) in [BZ_1979]
 double Atom::GetJdd(){
-   return this->zindoF0dd - 2.0*(this->zindoF2dd - this->zindoF4dd)/63.0;
+   return this->zindoF0dd - 2.0*(this->zindoF2dd + this->zindoF4dd)/63.0;
 }
 
 // Eq. (13) in [BZ_1979]
@@ -313,6 +363,123 @@ double Atom::GetIndoG1(){
 }
 
 
+// Table 1 in ref. [RZ_1976], Table 1 in [AEZ_1986], or Table 1 in [GD_1972]
+double Atom::GetZindoF0ss(){
+   return this->zindoF0ss;
+}
+
+// Table 1 in [AEZ_1986]
+double Atom::GetZindoF0sd(){
+   return this->zindoF0sd;
+}
+
+
+// Table 1 in [AEZ_1986]
+double Atom::GetZindoF0dd(){
+   return this->zindoF0dd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoG1sp(){
+   return this->zindoG1sp;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoF2pp(){
+   return this->zindoF2pp;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoG2sd(){
+   return this->zindoG2sd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoG1pd(){
+   return this->zindoG1pd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoF2pd(){
+   return this->zindoF2pd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoG3pd(){
+   return this->zindoG3pd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoF2dd(){
+   return this->zindoF2dd;
+}
+
+// Table 3 in ref. [BZ_1979]
+double Atom::GetZindoF4dd(){
+   return this->zindoF4dd;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF0ssLower(){
+   return this->zindoF0ss;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF0sdLower(){
+   return this->zindoF0sd;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF0ddLower(){
+   return this->zindoF0dd;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoG1spLower(){
+   return this->zindoG1sp/3.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF2ppLower(){
+   return this->zindoF2pp/25.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoG2sdLower(){
+   return this->zindoG2sd/5.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoG1pdLower(){
+   return this->zindoG1pd/15.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF2pdLower(){
+   return this->zindoF2pd/35.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoG3pdLower(){
+   return this->zindoG3pd/245.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF2ddLower(){
+   return this->zindoF2dd/49.0;
+}
+
+// Apendix in ref. [BZ_1979]
+double Atom::GetZindoF4ddLower(){
+   return this->zindoF4dd/441.0;
+}
+
+
+double Atom::GetCoreIntegral(OrbitalType orbital, bool isGuess, TheoryType theory){
+   return this->GetCoreIntegral(orbital, 0.0, isGuess, theory);
+}
+
+
 
 }
 #endif
@@ -321,3 +488,18 @@ double Atom::GetIndoG1(){
 
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

--- a/src/base/atoms/Catom.h
+++ b/src/base/atoms/Catom.h
@@ -21,6 +21,8 @@ Catom::Catom(double x, double y, double z) : Atom(x, y, z){
    this->valence.push_back(pz);
    this->valence.push_back(px);
    this->bondingParameter = -21.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterSZindo = -17.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterDZindo = 0.0;
    this->coreCharge = 4.0;
    this->imuAmuS = 14.051*Parameters::GetInstance()->GetEV2AU();
    this->imuAmuP = 5.572*Parameters::GetInstance()->GetEV2AU();
@@ -44,8 +46,8 @@ Catom::Catom(double x, double y, double z) : Atom(x, y, z){
    this->zindoG3pd = 0.0;                 
    this->zindoF2dd = 0.0;                 
    this->zindoF4dd = 0.0;                 
-   this->IonPotS = -19.84 * Parameters::GetInstance()->GetEV2AU();
-   this->IonPotP = -10.93 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotS = 19.84 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotP = 10.93 * Parameters::GetInstance()->GetEV2AU();
    this->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }
 

--- a/src/base/atoms/Hatom.h
+++ b/src/base/atoms/Hatom.h
@@ -19,6 +19,8 @@ Hatom::Hatom(double x, double y, double z) : Atom(x, y, z){
    this->atomType = H;
    this->valence.push_back(s);
    this->bondingParameter = -9.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterSZindo = -9.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterDZindo = 0.0;
    this->coreCharge = 1.0;
    this->imuAmuS = 7.176*Parameters::GetInstance()->GetEV2AU();
    this->imuAmuP = 0.0;
@@ -42,7 +44,7 @@ Hatom::Hatom(double x, double y, double z) : Atom(x, y, z){
    this->zindoG3pd = 0.0;                 
    this->zindoF2dd = 0.0;                 
    this->zindoF4dd = 0.0;                 
-   this->IonPotS = -13.06 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotS = 13.06 * Parameters::GetInstance()->GetEV2AU();
    this->IonPotP = 0.0 * Parameters::GetInstance()->GetEV2AU();
    this->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }

--- a/src/base/atoms/Liatom.h
+++ b/src/base/atoms/Liatom.h
@@ -22,6 +22,8 @@ Liatom::Liatom(double x, double y, double z) : Atom(x, y, z){
    this->valence.push_back(pz);
    this->valence.push_back(px);
    this->bondingParameter = -9.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterSZindo = 0.0;
+   this->bondingParameterDZindo = 0.0;
    this->coreCharge = 1.0;
    this->imuAmuS = 3.106*Parameters::GetInstance()->GetEV2AU();
    this->imuAmuP = 1.258*Parameters::GetInstance()->GetEV2AU();
@@ -45,8 +47,8 @@ Liatom::Liatom(double x, double y, double z) : Atom(x, y, z){
    this->zindoG3pd = 0.0;             
    this->zindoF2dd = 0.0;             
    this->zindoF4dd = 0.0;       
-   this->IonPotS = -5.39 * Parameters::GetInstance()->GetEV2AU();
-   this->IonPotP = -3.54 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotS = 5.39 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotP = 3.54 * Parameters::GetInstance()->GetEV2AU();
    this->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }
 

--- a/src/base/atoms/Satom.h
+++ b/src/base/atoms/Satom.h
@@ -29,6 +29,8 @@ Satom::Satom(double x, double y, double z) : Atom(x, y, z){
       this->valence.push_back(dxxyy);
    }
    this->bondingParameter = -18.150*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterSZindo = -14.0*Parameters::GetInstance()->GetEV2AU();
+   this->bondingParameterDZindo =   4.0*Parameters::GetInstance()->GetEV2AU();
    this->coreCharge = 6.0;
    this->imuAmuS = 17.650*Parameters::GetInstance()->GetEV2AU();
    this->imuAmuP = 6.989*Parameters::GetInstance()->GetEV2AU();
@@ -52,9 +54,9 @@ Satom::Satom(double x, double y, double z) : Atom(x, y, z){
    this->zindoG3pd = 20587*Parameters::GetInstance()->GetKayser2AU();           
    this->zindoF2dd = 28411*Parameters::GetInstance()->GetKayser2AU();           
    this->zindoF4dd = 18529*Parameters::GetInstance()->GetKayser2AU();           
-   this->IonPotS = -21.11 * Parameters::GetInstance()->GetEV2AU();
-   this->IonPotP = -12.39 * Parameters::GetInstance()->GetEV2AU();
-   this->IonPotD = -4.11 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotS = 21.11 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotP = 12.39 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotD = 4.11 * Parameters::GetInstance()->GetEV2AU();
 }
 
 double Satom::GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory){

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -57,12 +57,9 @@ private:
                                      double** orbitalElectronPopulation, Molecule* molecule);
    void CalcOverlap(double** overlap, Molecule* molecule);
    void CalcRotatingMatrix(double** rotatingMatrix, Atom* atomA, Atom* atomB);
-   void CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB);
    void CalcFockMatrix(double** fockMatrix, Molecule* molecule, double** overlap, double** gammaAB,
                        double** orbitalElectronPopulation, double* atomicElectronPopulation,
                        bool isGuess);
-   double GetReducedOverlap(int na, int la, int m, int nb, int lb, double alpha, double beta);
-   double GetReducedOverlap(int na, int nb, double alpha, double beta);
    void RotateDiatmicOverlapToSpaceFrame(double** diatomicOverlap, double** rotatingMatrix);
    void SetOverlapElement(double** overlap, double** diatomicOverlap, Atom* atomA, Atom* atomB);
    double GetAuxiliaryA(int k, double rho);
@@ -86,17 +83,19 @@ protected:
    string messageStartSCF;
    string messageDoneSCF;
    vector<AtomType> enableAtomTypes;
+   double GetReducedOverlap(int na, int la, int m, int nb, int lb, double alpha, double beta);
+   double GetReducedOverlap(int na, int nb, double alpha, double beta);
    virtual void CalcGammaAB(double** gammaAB, Molecule* molecule);
    virtual void SetMessages();
    virtual void SetEnableAtomTypes();
-   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
+   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, 
                              int mu, Molecule* molecule, double** gammaAB,
                              double** orbitalElectronPopulation, double* atomicElectronPopulation,
                              bool isGuess);
    virtual double GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                int firstAOIndexA, int firstAOIndexB,
                                 int mu, int nu, Molecule* molecule, double** gammaAB, double** overlap,
                                 double** orbitalElectronPopulation, bool isGuess);
+   virtual void CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB);
    TheoryType theory;
 public:
    Cndo2();
@@ -453,7 +452,6 @@ void Cndo2::CalcFockMatrix(double** fockMatrix, Molecule* molecule, double** ove
                   // diagonal part
                   fockMatrix[mu][mu] = this->GetFockDiagElement(atomA, 
                                                                 A, 
-                                                                firstAOIndexA, 
                                                                 mu, 
                                                                 molecule, 
                                                                 gammaAB,
@@ -467,8 +465,6 @@ void Cndo2::CalcFockMatrix(double** fockMatrix, Molecule* molecule, double** ove
                                                                    atomB,
                                                                    A, 
                                                                    B, 
-                                                                   firstAOIndexA, 
-                                                                   firstAOIndexB,
                                                                    mu, 
                                                                    nu, 
                                                                    molecule, 
@@ -499,11 +495,12 @@ void Cndo2::CalcFockMatrix(double** fockMatrix, Molecule* molecule, double** ove
 
 }
 
-double Cndo2::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, int mu, 
+double Cndo2::GetFockDiagElement(Atom* atomA, int atomAIndex, int mu, 
                                  Molecule* molecule, double** gammaAB,
                                  double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                  bool isGuess){
    double value;
+   int firstAOIndexA = atomA->GetFirstAOIndex();
    value = -1.0 * atomA->GetImuAmu(atomA->GetValence()[mu-firstAOIndexA]);
    if(!isGuess){
       double temp = atomicElectronPopulation[atomAIndex] - atomA->GetCoreCharge() 
@@ -525,7 +522,6 @@ double Cndo2::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA,
 }
 
 double Cndo2::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                    int firstAOIndexA, int firstAOIndexB,
                                     int mu, int nu, Molecule* molecule, double** gammaAB, double** overlap,
                                     double** orbitalElectronPopulation, bool isGuess){
    double value;

--- a/src/indo/Indo.h
+++ b/src/indo/Indo.h
@@ -23,15 +23,14 @@ private:
    double GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom); // Indo Exchange Interaction, (3.87) - (3.91) in J. A. Pople book.
 protected:
    virtual void SetMessages();
-   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
+   virtual void SetEnableAtomTypes();
+   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, 
                                      int mu, Molecule* molecule, double** gammaAB,
                                      double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                      bool isGuess);
    virtual double GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                int firstAOIndexA, int firstAOIndexB,
                                 int mu, int nu, Molecule* molecule, double** gammaAB, double** overelap,
                                 double** orbitalElectronPopulation, bool isGuess);
-   virtual void SetEnableAtomTypes();
 
 public:
    Indo();
@@ -77,11 +76,12 @@ void Indo::SetEnableAtomTypes(){
    this->enableAtomTypes.push_back(F);
 }
 
-double Indo::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, int mu, 
+double Indo::GetFockDiagElement(Atom* atomA, int atomAIndex, int mu, 
                                  Molecule* molecule, double** gammaAB,
                                  double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                  bool isGuess){
    double value;
+   int firstAOIndexA = atomA->GetFirstAOIndex();
    value = atomA->GetCoreIntegral(atomA->GetValence()[mu-firstAOIndexA], 
                                      gammaAB[atomAIndex][atomAIndex], 
                                      isGuess, this->theory);
@@ -102,11 +102,11 @@ double Indo::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA,
       value += temp;
    
       temp = 0.0;
-      for(int BB=0; BB<molecule->GetAtomVect()->size(); BB++){
-         if(BB != atomAIndex){
-            Atom* atomBB = (*molecule->GetAtomVect())[BB];
-            temp += ( atomicElectronPopulation[BB] - atomBB->GetCoreCharge()  )
-                     *gammaAB[atomAIndex][BB];
+      for(int B=0; B<molecule->GetAtomVect()->size(); B++){
+         if(B != atomAIndex){
+            Atom* atomB = (*molecule->GetAtomVect())[B];
+            temp += ( atomicElectronPopulation[B] - atomB->GetCoreCharge()  )
+                     *gammaAB[atomAIndex][B];
          }
       }
       value += temp;
@@ -116,7 +116,6 @@ double Indo::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA,
 }
 
 double Indo::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                    int firstAOIndexA, int firstAOIndexB,
                                     int mu, int nu, Molecule* molecule, double** gammaAB, double** overlap,
                                     double** orbitalElectronPopulation, bool isGuess){
    double value;
@@ -133,8 +132,8 @@ double Indo::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int
       double coulomb = 0.0;
       double exchange = 0.0;
       if(atomAIndex == atomBIndex){
-         OrbitalType orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
-         OrbitalType orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
+         OrbitalType orbitalMu = atomA->GetValence()[mu-atomA->GetFirstAOIndex()];
+         OrbitalType orbitalNu = atomA->GetValence()[nu-atomA->GetFirstAOIndex()];
          coulomb  = this->GetCoulombInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex], atomA); 
          exchange = this->GetExchangeInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex], atomA); 
          value = (1.5*exchange - 0.5*coulomb)*orbitalElectronPopulation[mu][nu];

--- a/src/input.in
+++ b/src/input.in
@@ -5,18 +5,18 @@ SCF_END
 
 THEORY
    //cndo/2
-   indo
-   //zindo/s
+   //indo
+   zindo/s
 THEORY_END
 
 //metane
-GEOMETRY
-   C  -0.37687006  0.95490165  0.00000000
-   H  -0.02021563  -0.05390835  0.00000000
-   H  -0.02019722  1.45929984  0.87365150
-   H  -0.02019722  1.45929984  -0.87365150
-   H  -1.44687006  0.95491483  0.00000000
-GEOMETRY_END
+//GEOMETRY
+//   C  -0.37687006  0.95490165  0.00000000
+//   H  -0.02021563  -0.05390835  0.00000000
+//   H  -0.02019722  1.45929984  0.87365150
+//   H  -0.02019722  1.45929984  -0.87365150
+//   H  -1.44687006  0.95491483  0.00000000
+//GEOMETRY_END
 
 
 // c2
@@ -38,11 +38,11 @@ GEOMETRY_END
 //GEOMETRY_END
 
 // sh2
-//GEOMETRY
-//     S       -0.559299    0.471698    0.000000
-//     H        0.750701    0.471698    0.000000
-//     H       -0.996586    1.706558    0.000000
-//GEOMETRY_END
+GEOMETRY
+     S       -0.559299    0.471698    0.000000
+     H        0.750701    0.471698    0.000000
+     H       -0.996586    1.706558    0.000000
+GEOMETRY_END
 
 // benzene
 //GEOMETRY

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -19,20 +19,27 @@ namespace MolDS_zindo{
  */
 class ZindoS : public MolDS_cndo::Cndo2{
 private:
-   double GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom); // Apendix in [BZ_1979]
-   double GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom); // Apendix in [BZ_1979]
+   double GetCoulombInt(OrbitalType orbital1, 
+                        OrbitalType orbital2, 
+                        Atom* atom); // Apendix in [BZ_1979]
+   double GetExchangeInt(OrbitalType orbital1, 
+                         OrbitalType orbital2, 
+                         Atom* atom); // Apendix in [BZ_1979]
+   double GetNishimotoMatagaTwoEleInt(Atom* atomA, OrbitalType orbitalA, 
+                                           Atom* atomB, OrbitalType orbitalB); // ref. [MN_1957] and (5a) in [AEZ_1986]
+   string errorMessageNishimotoMataga;
 protected:
    virtual void CalcGammaAB(double** gammaAB, Molecule* molecule);
    virtual void SetMessages();
-   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
+   virtual void SetEnableAtomTypes();
+   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, 
                                      int mu, Molecule* molecule, double** gammaAB,
                                      double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                      bool isGuess);
    virtual double GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                        int firstAOIndexA, int firstAOIndexB,
                                         int mu, int nu, Molecule* molecule, double** gammaAB, double** overelap,
                                         double** orbitalElectronPopulation, bool isGuess);
-   virtual void SetEnableAtomTypes();
+   virtual void CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB);
 public:
    ZindoS();
    ~ZindoS();
@@ -60,6 +67,7 @@ void ZindoS::SetMessages(){
       = "Error in zindo::ZindoS::CheckEnableAtomType: Non available atom is contained.\n";
    this->errorMessageCoulombInt = "Error in base_zindo::ZindoS::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_zindo::ZindoS::GetExchangeInt: Invalid orbitalType.\n";
+   this->errorMessageNishimotoMataga = "Error in base_zindo::ZindoS::GetNishimotoMatagaTwoEleInt: Invalid orbitalType.\n";
    this->messageSCFMetConvergence = "\n\n\n\t\tZINDO/S-SCF met convergence criterion(^^b\n\n\n";
    this->messageStartSCF = "**********  START: ZINDO/S-SCF  **********\n";
    this->messageDoneSCF = "**********  DONE: ZINDO/S-SCF  **********\n\n\n";
@@ -74,14 +82,14 @@ void ZindoS::SetEnableAtomTypes(){
    this->enableAtomTypes.push_back(S);
 }
 
-double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, int mu, 
+double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int mu, 
                                  Molecule* molecule, double** gammaAB,
                                  double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                  bool isGuess){
-   double value;
+   double value=0.0;
+   int firstAOIndexA = atomA->GetFirstAOIndex();
    value = atomA->GetCoreIntegral(atomA->GetValence()[mu-firstAOIndexA], 
-                                     gammaAB[atomAIndex][atomAIndex], 
-                                     isGuess, this->theory);
+                                  isGuess, this->theory);
    if(!isGuess){
       double temp = 0.0;
       double coulomb = 0.0;
@@ -90,19 +98,25 @@ double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA
       OrbitalType orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
       for(int v=0; v<atomA->GetValence().size(); v++){
          OrbitalType orbitalLam = atomA->GetValence()[v];
-         coulomb  = this->GetCoulombInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
-         exchange = this->GetExchangeInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
-         lammda = firstAOIndexA + v;
+         coulomb  = this->GetCoulombInt(orbitalMu, orbitalLam, atomA);
+         exchange = this->GetExchangeInt(orbitalMu, orbitalLam, atomA);
+         lammda = v + firstAOIndexA;
          temp += orbitalElectronPopulation[lammda][lammda]*(coulomb - 0.5*exchange);
       }
       value += temp;
    
       temp = 0.0;
-      for(int BB=0; BB<molecule->GetAtomVect()->size(); BB++){
-         if(BB != atomAIndex){
-            Atom* atomBB = (*molecule->GetAtomVect())[BB];
-            temp += ( atomicElectronPopulation[BB] - atomBB->GetCoreCharge()  )
-                     *gammaAB[atomAIndex][BB];
+      for(int B=0; B<molecule->GetAtomVect()->size(); B++){
+         if(B != atomAIndex){
+            Atom* atomB = (*molecule->GetAtomVect())[B];
+            for(int i=0; i<atomB->GetValence().size(); i++){
+               int sigma = i + atomB->GetFirstAOIndex();
+               OrbitalType orbitalSigma = atomB->GetValence()[i];
+               temp += orbitalElectronPopulation[sigma][sigma]
+                      *this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalSigma);
+            }
+            temp -= atomB->GetCoreCharge() 
+                   *this->GetNishimotoMatagaTwoEleInt(atomA, s, atomB, s);
          }
       }
       value += temp;
@@ -112,15 +126,14 @@ double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA
 }
 
 double ZindoS::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
-                                    int firstAOIndexA, int firstAOIndexB,
                                     int mu, int nu, Molecule* molecule, double** gammaAB, double** overlap,
                                     double** orbitalElectronPopulation, bool isGuess){
-   double value;
-   double K = 1.0;
-   if(m <= atomA->GetValenceShellType() || m <= atomB->GetValenceShellType()){
-      K = 0.75;
-   }
-   double bondParameter = 0.5*K*(atomA->GetBondingParameter() + atomB->GetBondingParameter()); 
+   double value = 0.0;
+   OrbitalType orbitalMu = atomA->GetValence()[mu-atomA->GetFirstAOIndex()];
+   OrbitalType orbitalNu = atomB->GetValence()[nu-atomB->GetFirstAOIndex()];
+
+   double bondParameter = 0.5*(atomA->GetBondingParameter(this->theory, orbitalMu) 
+                              +atomB->GetBondingParameter(this->theory, orbitalNu)); 
 
    if(isGuess){
       value = bondParameter*overlap[mu][nu];
@@ -129,18 +142,17 @@ double ZindoS::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, i
       double coulomb = 0.0;
       double exchange = 0.0;
       if(atomAIndex == atomBIndex){
-         OrbitalType orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
-         OrbitalType orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
-         coulomb  = this->GetCoulombInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex], atomA); 
-         exchange = this->GetExchangeInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex], atomA); 
+         coulomb  = this->GetCoulombInt(orbitalMu, orbitalNu, atomA); 
+         exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA); 
          value = (1.5*exchange - 0.5*coulomb)*orbitalElectronPopulation[mu][nu];
       }
       else{
          value = bondParameter*overlap[mu][nu];
-         value -= 0.5*orbitalElectronPopulation[mu][nu]*gammaAB[atomAIndex][atomBIndex];
+         value -= 0.5*orbitalElectronPopulation[mu][nu]
+                  *this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalNu);
       }
    }
-
+   
    return value;
 }
 
@@ -150,29 +162,131 @@ void ZindoS::CalcGammaAB(double** gammaAB, Molecule* molecule){
 
 // Apendix in [BZ_1972]
 // ZINDO Coulomb Interaction
-double ZindoS::GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+double ZindoS::GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, Atom* atom){
 
    double value=0.0;
-
-   // ToDo: Coulomb interaction
-   /*
+   
    if( orbital1 == s && orbital2 == s){ 
-      value = gamma;
+      value = atom->GetZindoF0ssLower();
    }   
    else if( orbital1 == s && ( orbital2 == px || orbital2 == py || orbital2 == pz )){ 
-      value = gamma;
+      value = atom->GetZindoF0ssLower();
    }   
-   else if( (orbital1 == px || orbital1 == py || orbital1 == pz ) && orbital2 == s){ 
-      value = gamma;
+   else if( orbital2 == s && ( orbital1 == px || orbital1 == py || orbital1 == pz )){ 
+      value = atom->GetZindoF0ssLower();
    }   
    else if( (orbital1 == orbital2) && ( orbital1 == px || orbital1 == py || orbital1 == pz )){ 
-      value = gamma + 4.0*atom->GetIndoF2()/25.0;
+      value = atom->GetZindoF0ssLower()
+             +atom->GetZindoF2ppLower()*4.0;
    }   
    else if( (orbital1 != orbital2) 
          && ( orbital1 == px || orbital1 == py || orbital1 == pz )
          && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = gamma - 2.0*atom->GetIndoF2()/25.0;
+      value = atom->GetZindoF0ssLower()
+             -atom->GetZindoF2ppLower()*2.0;
+   }   
+   else if( orbital1 == s && ( orbital2 == dxy || 
+                               orbital2 == dyz || 
+                               orbital2 == dzz || 
+                               orbital2 == dzx || 
+                               orbital2 == dxxyy )){ 
+      value = atom->GetZindoF0sdLower();
    }   
+   else if( orbital2 == s && ( orbital1 == dxy || 
+                               orbital1 == dyz || 
+                               orbital1 == dzz || 
+                               orbital1 == dzx || 
+                               orbital1 == dxxyy )){ 
+      value = atom->GetZindoF0sdLower();
+   }
+   else if( orbital1 == dzz && (orbital2 == px || orbital2==py) ){
+      value = atom->GetZindoF0sdLower()
+             -atom->GetZindoF2pdLower()*2.0;
+   }
+   else if( orbital2 == dzz && (orbital1 == px || orbital1==py) ){
+      value = atom->GetZindoF0sdLower()
+             -atom->GetZindoF2pdLower()*2.0;
+   }
+   else if( (orbital1 == dzz && orbital2 == pz) ||
+            (orbital2 == dzz && orbital1 == pz) ){
+      value = atom->GetZindoF0sdLower()
+             +atom->GetZindoF2pdLower()*4.0;
+   }
+   else if( (orbital1 == orbital2) && ( orbital1 == dxy || 
+                                        orbital1 == dyz || 
+                                        orbital1 == dzz || 
+                                        orbital1 == dzx || 
+                                        orbital1 == dxxyy )){ 
+      value = atom->GetZindoF0ddLower()
+             +atom->GetZindoF2ddLower()*4.0
+             +atom->GetZindoF4ddLower()*36.0;
+   }
+   else if( (orbital1 == dxxyy && orbital2 == px) ||
+            (orbital2 == dxxyy && orbital1 == px) || 
+            (orbital1 == dxxyy && orbital2 == py) ||
+            (orbital2 == dxxyy && orbital1 == py) || 
+            (orbital1 == dxy && orbital2 == px) ||
+            (orbital2 == dxy && orbital1 == px) || 
+            (orbital1 == dxy && orbital2 == py) ||
+            (orbital2 == dxy && orbital1 == py) ||
+            (orbital1 == dzx && orbital2 == px) ||
+            (orbital2 == dzx && orbital1 == px) || 
+            (orbital1 == dzx && orbital2 == pz) ||
+            (orbital2 == dzx && orbital1 == pz) || 
+            (orbital1 == dyz && orbital2 == py) ||
+            (orbital2 == dyz && orbital1 == py) || 
+            (orbital1 == dyz && orbital2 == pz) ||
+            (orbital2 == dyz && orbital1 == pz) ){
+      value = atom->GetZindoF0sdLower()
+             +atom->GetZindoF2pdLower()*2.0;
+   }
+   else if( (orbital1 == dxxyy && orbital2 == pz) ||
+            (orbital2 == dxxyy && orbital1 == pz) || 
+            (orbital1 == dxy && orbital2 == pz) ||
+            (orbital2 == dxy && orbital1 == pz) ||
+            (orbital1 == dzx && orbital2 == py) ||
+            (orbital2 == dzx && orbital1 == py) ||
+            (orbital1 == dyz && orbital2 == px) ||
+            (orbital2 == dyz && orbital1 == px)  ){
+      value = atom->GetZindoF0sdLower()
+             -atom->GetZindoF2pdLower()*4.0;
+   }
+   else if( (orbital1 == dxxyy && orbital2 == dzz) ||
+            (orbital2 == dxxyy && orbital1 == dzz) || 
+            (orbital1 == dxy && orbital2 == dzz) ||
+            (orbital2 == dxy && orbital1 == dzz)  ){
+      value = atom->GetZindoF0ddLower()
+             -atom->GetZindoF2ddLower()*4.0
+             +atom->GetZindoF4ddLower()*6.0;
+   }
+   else if( (orbital1 == dxy && orbital2 == dxxyy) ||
+            (orbital2 == dxy && orbital1 == dxxyy)  ){
+      value = atom->GetZindoF0ddLower()
+             +atom->GetZindoF2ddLower()*4.0
+             -atom->GetZindoF4ddLower()*34.0;
+   }
+   else if( (orbital1 == dzx && orbital2 == dzz) ||
+            (orbital2 == dzx && orbital1 == dzz) || 
+            (orbital1 == dyz && orbital2 == dzz) ||
+            (orbital2 == dyz && orbital1 == dzz)  ){
+      value = atom->GetZindoF0ddLower()
+             +atom->GetZindoF2ddLower()*2.0
+             -atom->GetZindoF4ddLower()*24.0;
+   }
+   else if( (orbital1 == dzx && orbital2 == dxxyy) ||
+            (orbital2 == dzx && orbital1 == dxxyy) || 
+            (orbital1 == dzx && orbital2 == dxy) || 
+            (orbital2 == dzx && orbital1 == dxy) || 
+            (orbital1 == dyz && orbital2 == dxxyy) || 
+            (orbital2 == dyz && orbital1 == dxxyy) || 
+            (orbital1 == dyz && orbital2 == dxy) || 
+            (orbital2 == dyz && orbital1 == dxy) || 
+            (orbital1 == dyz && orbital2 == dzx) || 
+            (orbital2 == dyz && orbital1 == dzx) ){
+      value = atom->GetZindoF0ddLower()
+             -atom->GetZindoF2ddLower()*2.0
+             -atom->GetZindoF4ddLower()*4.0;
+   }
    else{
       stringstream ss;
       ss << this->errorMessageCoulombInt;
@@ -181,33 +295,117 @@ double ZindoS::GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double
       ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
       throw MolDSException(ss.str());
    }   
-   */
+   
    return value;
 
 }
 
 // Apendix in [BZ_1972]
 // ZINDO Exchange Interaction
-double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, Atom* atom){
 
    double value=0.0;
 
-   // ToDo: Exchange interaction
-   /*
    if( orbital1 == orbital2){
-      value = this->GetCoulombInt(orbital1, orbital2, gamma, atom);
+      value = this->GetCoulombInt(orbital1, orbital2, atom);
+   }   
+   else if( orbital1 == s && (orbital2 == px || orbital2 == py || orbital2 == pz ) ){
+      value = atom->GetZindoG1spLower();
+   }   
+   else if( orbital2 == s && (orbital1 == px || orbital1 == py || orbital1 == pz ) ){
+      value = atom->GetZindoG1spLower();
    }   
-   else if( (orbital1 == s) && (orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = atom->GetIndoG1()/3.0;
+   else if( (orbital1 == s) && (orbital2 == dxy || 
+                                orbital2 == dyz || 
+                                orbital2 == dzz || 
+                                orbital2 == dzx || 
+                                orbital2 == dxxyy ) ){
+      value = atom->GetZindoG2sdLower();
    }   
-   else if( (orbital1 == px || orbital1 == py || orbital1 == pz) && orbital2 == s  ){  
-      value = atom->GetIndoG1()/3.0;
+   else if( (orbital2 == s) && (orbital1 == dxy || 
+                                orbital1 == dyz || 
+                                orbital1 == dzz || 
+                                orbital1 == dzx || 
+                                orbital1 == dxxyy ) ){
+      value = atom->GetZindoG2sdLower();
    }   
    else if( (orbital1 != orbital2) 
          && ( orbital1 == px || orbital1 == py || orbital1 == pz )
          && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = 3.0*atom->GetIndoF2()/25.0;
-   }   
+      value = atom->GetZindoF2ppLower()*3.0;
+   }
+   else if( (orbital1 == px && orbital2 == dzz) ||
+            (orbital2 == px && orbital1 == dzz) ||
+            (orbital1 == py && orbital2 == dzz) ||
+            (orbital2 == py && orbital1 == dzz) ){
+      value = atom->GetZindoG1pdLower()
+             +atom->GetZindoG3pdLower()*18.0;
+   }
+   else if( (orbital1 == px && orbital2 == dxxyy) ||
+            (orbital2 == px && orbital1 == dxxyy) ||
+            (orbital1 == px && orbital2 == dxy) ||
+            (orbital2 == px && orbital1 == dxy) ||
+            (orbital1 == px && orbital2 == dzx) ||
+            (orbital2 == px && orbital1 == dzx) ||
+            (orbital1 == py && orbital2 == dxxyy) ||
+            (orbital2 == py && orbital1 == dxxyy) ||
+            (orbital1 == py && orbital2 == dxy) ||
+            (orbital2 == py && orbital1 == dxy) ||
+            (orbital1 == py && orbital2 == dyz) ||
+            (orbital2 == py && orbital1 == dyz) ||
+            (orbital1 == pz && orbital2 == dzx) ||
+            (orbital2 == pz && orbital1 == dzx) ||
+            (orbital1 == pz && orbital2 == dyz) ||
+            (orbital2 == pz && orbital1 == dyz) ){
+      value = atom->GetZindoG1pdLower()*3.0
+             +atom->GetZindoG3pdLower()*24.0;
+   }
+   else if( (orbital1 == px && orbital2 == dyz) ||
+            (orbital2 == px && orbital1 == dyz) ||
+            (orbital1 == py && orbital2 == dzx) ||
+            (orbital2 == py && orbital1 == dzx) ||
+            (orbital1 == pz && orbital2 == dxxyy) ||
+            (orbital2 == pz && orbital1 == dxxyy) ||
+            (orbital1 == pz && orbital2 == dxy) ||
+            (orbital2 == pz && orbital1 == dxy) ){
+      value = atom->GetZindoG3pdLower()*15.0;
+   }
+   else if( (orbital1 == pz && orbital2 == dzz) ||
+            (orbital2 == pz && orbital1 == dzz) ){
+      value = atom->GetZindoG1pdLower()*4.0
+             +atom->GetZindoG3pdLower()*27.0;
+   }
+   else if( (orbital1 == dzz && orbital2 == dxxyy) ||
+            (orbital2 == dzz && orbital1 == dxxyy) ||
+            (orbital1 == dzz && orbital2 == dxy) ||
+            (orbital2 == dzz && orbital1 == dxy) ){
+      value = atom->GetZindoF2ddLower()*4.0
+             +atom->GetZindoF4ddLower()*15.0;
+   }
+   else if( (orbital1 == dzz && orbital2 == dzx) ||
+            (orbital2 == dzz && orbital1 == dzx) ||
+            (orbital1 == dzz && orbital2 == dyz) ||
+            (orbital2 == dzz && orbital1 == dyz) ){
+      value = atom->GetZindoF2ddLower()
+             +atom->GetZindoF4ddLower()*30.0;
+   }
+   else if( (orbital1 == dxxyy && orbital2 == dxy) ||
+            (orbital2 == dxxyy && orbital1 == dxy) ){
+      value = atom->GetZindoF4ddLower()*35.0;
+   }
+   else if( (orbital1 == dxxyy && orbital2 == dzx) ||
+            (orbital2 == dxxyy && orbital1 == dzx) ||
+            (orbital1 == dxxyy && orbital2 == dyz) ||
+            (orbital2 == dxxyy && orbital1 == dyz) ||
+            (orbital1 == dxy && orbital2 == dzx) ||
+            (orbital2 == dxy && orbital1 == dzx) ||
+            (orbital1 == dxy && orbital2 == dyz) ||
+            (orbital2 == dxy && orbital1 == dyz) ||
+            (orbital1 == dzx && orbital2 == dyz) ||
+            (orbital2 == dzx && orbital1 == dyz) ){
+      value = atom->GetZindoF2ddLower()*3.0
+             +atom->GetZindoF4ddLower()*20.0;
+   }
    else{
       stringstream ss;
       ss << this->errorMessageExchangeInt;
@@ -216,11 +414,92 @@ double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double
       ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
       throw MolDSException(ss.str());
    }   
-   */
+   
 
    return value;
 }
 
+// ref. [MN_1957] and (5a) in [AEZ_1986]
+double ZindoS::GetNishimotoMatagaTwoEleInt(Atom* atomA, OrbitalType orbitalA, 
+                                           Atom* atomB, OrbitalType orbitalB){
+
+   double r =sqrt(
+                   pow(atomA->GetXyz()[0] - atomB->GetXyz()[0], 2.0)
+                  +pow(atomA->GetXyz()[1] - atomB->GetXyz()[1], 2.0)
+                  +pow(atomA->GetXyz()[2] - atomB->GetXyz()[2], 2.0) );
+  
+   double gammaAA;
+   if(orbitalA == s || 
+      orbitalA == px ||
+      orbitalA == py ||
+      orbitalA == pz ){
+      gammaAA = atomA->GetZindoF0ss();
+   }
+   else if(orbitalA == dxy ||
+           orbitalA == dyz ||
+           orbitalA == dzz ||
+           orbitalA == dzx ||
+           orbitalA == dxxyy ){
+      gammaAA = atomA->GetZindoF0dd();
+   }
+   else{
+      stringstream ss;
+      ss << this->errorMessageNishimotoMataga;
+      ss << this->errorMessageAtomType << AtomTypeStr(atomA->GetAtomType()) << "\n";
+      ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalA) << "\n";
+      throw MolDSException(ss.str());
+   }   
+
+   double gammaBB;
+   if(orbitalB == s || 
+      orbitalB == px ||
+      orbitalB == py ||
+      orbitalB == pz ){
+      gammaBB = atomB->GetZindoF0ss();
+   }
+   else if(orbitalB == dxy ||
+           orbitalB == dyz ||
+           orbitalB == dzz ||
+           orbitalB == dzx ||
+           orbitalB == dxxyy ){
+      gammaBB = atomB->GetZindoF0dd();
+   }
+   else{
+      stringstream ss;
+      ss << this->errorMessageNishimotoMataga;
+      ss << this->errorMessageAtomType << AtomTypeStr(atomB->GetAtomType()) << "\n";
+      ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalB) << "\n";
+      throw MolDSException(ss.str());
+   }  
+
+   return 1.2/( r+2.4/(gammaAA+gammaBB) );
+
+}
+
+
+void ZindoS::CalcDiatomicOverlapInDiatomicFrame(double** diatomicOverlap, Atom* atomA, Atom* atomB){
+
+   MolDS_cndo::Cndo2::CalcDiatomicOverlapInDiatomicFrame(diatomicOverlap, atomA, atomB);
+
+   // see (4f) in [AEZ_1986]
+   diatomicOverlap[pz][pz] *= 1.267;
+   diatomicOverlap[py][py] *= 0.585;
+   diatomicOverlap[px][px] *= 0.585;
+
+   /*
+   for(int i=0;i<OrbitalType_end;i++){
+      for(int j=0;j<OrbitalType_end;j++){
+         printf("diatomicOverlap[%d][%d]=%lf\n",i,j,diatomicOverlap[i][j]);
+      }
+   }
+   */
+
+
+}
+
+
+
+
 }
 #endif