提交: c5ef493f4cc9b9ada4f1537c1ba51e75428db726 - MolDS (git)

Commit MetaInfo

• 修订版: c5ef493f4cc9b9ada4f1537c1ba51e75428db726 (tree)
• 时间: 2010-12-09 12:38:53
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Core and Exchnge integral for INDO and ZINDO are moved from Atom class to INDO and ZINDOS class, respectively.

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

更改概述

• modified: src/MolDS.cpp
• modified: src/base/atoms/Atom.h
• modified: src/cndo/Cndo2.h
• modified: src/indo/Indo.h
• modified: src/input.in
• modified: src/zindo/ZindoS.h

差异

--- a/src/MolDS.cpp
+++ b/src/MolDS.cpp
@@ -23,6 +23,7 @@
 #include"base/Parameters.h"
 #include"cndo/Cndo2.h"
 #include"indo/Indo.h"
+#include"zindo/ZindoS.h"
 #include"mkl_wrapper/LapackWrapper.h"
 
 
@@ -72,6 +73,10 @@ int main(){
 
    // ZINDO/S
    else if(Parameters::GetInstance()->GetCurrentTheory() == ZINDOS){
+      MolDS_zindo::ZindoS* zindoS = new MolDS_zindo::ZindoS();
+      //zindoS->SetMolecule(molecule);
+      //zindoS->DoesSCF();
+      delete zindoS;
    }
 
    /*** test lapack ***

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -17,8 +17,6 @@ private:
    void SetMessages();
    string errorMessageImuAmu;
    string errorMessageOrbitalExponent;
-   string errorMessageIndoCoulombInt;
-   string errorMessageIndoExchangeInt;
    string errorMessageShellType;
    string errorMessageEffectivPrincipalQuantumNumber;
    string errorMessageZindoCoreIntegral;
@@ -83,9 +81,9 @@ public:
    int GetNumberValenceElectrons();
    double GetImuAmu(OrbitalType orbitalType);  // return 0.5*(I_mu + A_mu)
    double GetOrbitalExponent(ShellType shellType, OrbitalType orbitalType);  // (1.73) in J. A. Pople book.
-   double GetIndoCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma); // (3.87) - (3.91) in J. A. Pople book.
-   double GetIndoExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma); // (3.87) - (3.91) 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 GetIndoF2();
+   double GetIndoG1();
 };
 
 Atom::Atom(){
@@ -111,8 +109,6 @@ Atom::~Atom(){
 void Atom::SetMessages(){
    this->errorMessageImuAmu = "Error in base_atoms::Atom::GetImuAmu: Invalid orbitalType.\n";
    this->errorMessageOrbitalExponent = "Error in base_atoms::Atom::GetOrbitalExponent: Invalid shelltype or orbitalType.\n";
-   this->errorMessageIndoCoulombInt = "Error in base_atoms::Atom::GetIndoCoulombInt: Invalid orbitalType.\n";
-   this->errorMessageIndoExchangeInt = "Error in base_atoms::Atom::GetIndoExchangeInt: Invalid orbitalType.\n";
    this->errorMessageIndoCoreIntegral = "Error in base_atoms::Atom::GetCoreIntegral: Invalid orbitalType for INDO.\n";
    this->errorMessageZindoSCoreIntegral = "Error in base_atoms::Atom::GetCoreIntegral: Invalid orbitalType for ZINDO/S.\n";
    this->errorMessageAtomType = "\tatom type = ";
@@ -239,69 +235,6 @@ double Atom::GetOrbitalExponent(ShellType shellType, OrbitalType orbitalType){
 }
 
 
-// (3.87) - (3.91) in J. A. Pople book.
-// Indo Coulomb Interaction
-double Atom::GetIndoCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma){
-
-   double value=0.0;
-   if( orbital1 == s && orbital2 == s){
-      value = gamma;
-   }
-   else if( orbital1 == s && ( orbital2 == px || orbital2 == py || orbital2 == pz )){
-      value = gamma;
-   }
-   else if( (orbital1 == px || orbital1 == py || orbital1 == pz ) && orbital2 == s){
-      value = gamma;
-   }
-   else if( (orbital1 == orbital2) && ( orbital1 == px || orbital1 == py || orbital1 == pz )){
-      value = gamma + 4.0*this->indoF2/25.0;
-   }
-   else if( (orbital1 != orbital2) 
-         && ( orbital1 == px || orbital1 == py || orbital1 == pz )
-         && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = gamma - 2.0*this->indoF2/25.0;
-   }
-   else{
-      cout << this->errorMessageIndoCoulombInt;
-      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << "\n";
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
-      exit(EXIT_FAILURE);
-   }
-
-   return value;
-}
-
-// (3.87) - (3.91) in J. A. Pople book.
-// Indo Exchange interaction
-double Atom::GetIndoExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma){
-   double value=0.0;
-
-   if( orbital1 == orbital2){
-      value = this->GetIndoCoulombInt(orbital1, orbital2, gamma);
-   }
-   else if( (orbital1 == s) && (orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = this->indoG1/3.0;
-   }
-   else if( (orbital1 == px || orbital1 == py || orbital1 == pz) && orbital2 == s  ){
-      value = this->indoG1/3.0;
-   }
-   else if( (orbital1 != orbital2) 
-         && ( orbital1 == px || orbital1 == py || orbital1 == pz )
-         && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
-      value = 3.0*this->indoF2/25.0;
-   }
-   else{
-      cout << this->errorMessageIndoExchangeInt;
-      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << "\n";
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
-      exit(EXIT_FAILURE);
-   }
-
-   return value;
-}
-
 // Part of Eq. (13) in [BZ_1979]
 double Atom::GetJss(){
    return this->zindoF0ss;
@@ -367,6 +300,13 @@ double Atom::GetZindoCoreIntegral(OrbitalType orbital, int l, int m, int n){
    return value;
 }
 
+double Atom::GetIndoF2(){
+   return this->indoF2;
+}
+
+double Atom::GetIndoG1(){
+   return this->indoG1;
+}
 
 
 

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -24,7 +24,6 @@ namespace MolDS_cndo{
  */
 class Cndo2{
 private:
-   string errorMessageAtomType;
    string messageEnergiesMOs;
    string messageEnergiesMOsTitle;
    string messageMullikenAtoms;
@@ -73,10 +72,14 @@ private:
    void CheckNumberValenceElectrons(Molecule* molecule);
 
 protected:
+   string errorMessageAtomType;
+   string errorMessageOrbitalType;
    string errorMessageSCFNotConverged;
    string errorMessageMoleculeNotSet;
    string errorMessageOddTotalValenceElectrions;
    string errorMessageNotEnebleAtomType;
+   string errorMessageCoulombInt;
+   string errorMessageExchangeInt;
    string messageSCFMetConvergence;
    string messageStartSCF;
    string messageDoneSCF;
@@ -156,8 +159,8 @@ void Cndo2::SetMessages(){
       = "Error in cndo::Cndo2::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in cndo::Cndo2::ChecEnableAtomType: Not enable atom is contained.\n";
-   this->errorMessageAtomType  
-      = "\tatom type = ";
+   this->errorMessageAtomType = "\tatom type = ";
+   this->errorMessageOrbitalType = "\torbital type = ";
    this->messageSCFMetConvergence = "\n\n\n\t\tCNDO/2-SCF met convergence criterion(^^b\n\n\n";
    this->messageStartSCF = "**********  START: CNDO/2-SCF  **********\n";
    this->messageDoneSCF = "**********  DONE: CNDO/2-SCF  **********\n\n\n";

--- a/src/indo/Indo.h
+++ b/src/indo/Indo.h
@@ -18,20 +18,24 @@ namespace MolDS_indo{
  *  Refferences for Indo are [PB_1970] and [PS_1966].
  */
 class Indo : public MolDS_cndo::Cndo2{
-public:
-   Indo();
-   ~Indo();
+private:
+   double GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom); // Indo Coulomb Interaction, (3.87) - (3.91) in J. A. Pople book.
+   double GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom); // Indo Exchange Interaction, (3.87) - (3.91) in J. A. Pople book.
 protected:
-   void SetMessages();
-   double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
-                             int mu, Molecule* molecule, double** gammaAB,
-                             double** orbitalElectronPopulation, double* atomicElectronPopulation,
-                             bool isGuess);
-   double GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int atomBIndex, 
+   virtual void SetMessages();
+   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
+                                     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);
-   void SetEnableAtomTypes();
+   virtual void SetEnableAtomTypes();
+
+public:
+   Indo();
+   ~Indo();
 };
 
 Indo::Indo() : MolDS_cndo::Cndo2(){
@@ -54,6 +58,8 @@ void Indo::SetMessages(){
       = "Error in indo::Indo::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in indo::Indo::ChecEnableAtomType: Not enable atom is contained.\n";
+   this->errorMessageCoulombInt = "Error in base_indo::Indo::GetCoulombInt: Invalid orbitalType.\n";
+   this->errorMessageExchangeInt = "Error in base_indo::Indo::GetExchangeInt: Invalid orbitalType.\n";
    this->messageSCFMetConvergence = "\n\n\n\t\tINDO-SCF met convergence criterion(^^b\n\n\n";
    this->messageStartSCF = "**********  START: INDO-SCF  **********\n";
    this->messageDoneSCF = "**********  DONE: INDO-SCF  **********\n\n\n";
@@ -88,8 +94,8 @@ double Indo::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  = atomA->GetIndoCoulombInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex]);
-         exchange = atomA->GetIndoExchangeInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex]);
+         coulomb  = this->GetCoulombInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
+         exchange = this->GetExchangeInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
          lammda = firstAOIndexA + v;
          temp += orbitalElectronPopulation[lammda][lammda]*(coulomb - 0.5*exchange);
       }
@@ -129,8 +135,8 @@ double Indo::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int
       if(atomAIndex == atomBIndex){
          OrbitalType orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
          OrbitalType orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
-         coulomb  = atomA->GetIndoCoulombInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex]); 
-         exchange = atomA->GetIndoExchangeInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex]); 
+         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];
       }
       else{
@@ -142,6 +148,71 @@ double Indo::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, int
    return value;
 }
 
+// (3.87) - (3.91) in J. A. Pople book.
+// Indo Coulomb Interaction
+double Indo::GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+
+   double value=0.0;
+   if( orbital1 == s && orbital2 == s){ 
+      value = gamma;
+   }   
+   else if( orbital1 == s && ( orbital2 == px || orbital2 == py || orbital2 == pz )){ 
+      value = gamma;
+   }   
+   else if( (orbital1 == px || orbital1 == py || orbital1 == pz ) && orbital2 == s){ 
+      value = gamma;
+   }   
+   else if( (orbital1 == orbital2) && ( orbital1 == px || orbital1 == py || orbital1 == pz )){ 
+      value = gamma + 4.0*atom->GetIndoF2()/25.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;
+   }   
+   else{
+      cout << this->errorMessageCoulombInt;
+      cout << this->errorMessageAtomType << AtomTypeStr(atom->GetAtomType()) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
+      exit(EXIT_FAILURE);
+   }   
+
+   return value;
+
+}
+
+// (3.87) - (3.91) in J. A. Pople book.
+// Indo Exchange Interaction
+double Indo::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+
+   double value=0.0;
+
+   if( orbital1 == orbital2){
+      value = this->GetCoulombInt(orbital1, orbital2, gamma, atom);
+   }   
+   else if( (orbital1 == s) && (orbital2 == px || orbital2 == py || orbital2 == pz ) ){
+      value = atom->GetIndoG1()/3.0;
+   }   
+   else if( (orbital1 == px || orbital1 == py || orbital1 == pz) && orbital2 == s  ){  
+      value = atom->GetIndoG1()/3.0;
+   }   
+   else if( (orbital1 != orbital2) 
+         && ( orbital1 == px || orbital1 == py || orbital1 == pz )
+         && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
+      value = 3.0*atom->GetIndoF2()/25.0;
+   }   
+   else{
+      cout << this->errorMessageExchangeInt;
+      cout << this->errorMessageAtomType << AtomTypeStr(atom->GetAtomType()) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
+      exit(EXIT_FAILURE);
+   }   
+
+   return value;
+
+}
 
 
 }

--- a/src/input.in
+++ b/src/input.in
@@ -45,29 +45,29 @@ THEORY_END
 //GEOMETRY_END
 
 // benzene
-//GEOMETRY
-//   C  -0.505391      0.552561      0.000000
-//   C  0.889769      0.552561      0.000000
-//   C  1.587307      1.760312      0.000000
-//   C  0.889653      2.968821     -0.001199
-//   C  -0.505172      2.968743     -0.001678
-//   C  -1.202773      1.760537     -0.000682
-//   H  -1.055150     -0.399756      0.000450
-//   H  1.439277     -0.399952      0.001315
-//   H  2.686987      1.760392      0.000634
-//   H  1.439853      3.920964     -0.001258
-//   H  -1.055294      3.921024     -0.002631
-//   H  -2.302377      1.760720     -0.000862
-//GEOMETRY_END
-
-// acetylene
 GEOMETRY
- C                     0.24933   0.0   0.0 
- C                     1.45053   0.0   0.0 
- H                    -0.82067   0.0   0.0 
- H                     2.52053   0.0   0.0 
+   C  -0.505391      0.552561      0.000000
+   C  0.889769      0.552561      0.000000
+   C  1.587307      1.760312      0.000000
+   C  0.889653      2.968821     -0.001199
+   C  -0.505172      2.968743     -0.001678
+   C  -1.202773      1.760537     -0.000682
+   H  -1.055150     -0.399756      0.000450
+   H  1.439277     -0.399952      0.001315
+   H  2.686987      1.760392      0.000634
+   H  1.439853      3.920964     -0.001258
+   H  -1.055294      3.921024     -0.002631
+   H  -2.302377      1.760720     -0.000862
 GEOMETRY_END
 
+// acetylene
+//GEOMETRY
+// C                     0.24933   0.0   0.0 
+// C                     1.45053   0.0   0.0 
+// H                    -0.82067   0.0   0.0 
+// H                     2.52053   0.0   0.0 
+//GEOMETRY_END
+
 // acetylene2
 //GEOMETRY
 // C           0.0, 0.0          0.24933 

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -1,5 +1,5 @@
-#ifndef INCLUDED_INDO
-#define INCLUDED_INDO
+#ifndef INCLUDED_ZINDOS
+#define INCLUDED_ZINDOS
 
 #include<stdio.h>
 #include<stdlib.h>
@@ -15,24 +15,27 @@ using namespace MolDS_base_atoms;
 namespace MolDS_zindo{
 
 /***
- *  Refferences for Indo are [PB_1970] and [PS_1966].
+ *  Refferences for 
  */
 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]
+protected:
+   virtual void CalcGammaAB(double** gammaAB, Molecule* molecule);
+   virtual void SetMessages();
+   virtual double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
+                                     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:
    ZindoS();
    ~ZindoS();
-protected:
-   void CalcGammaAB(double** gammaAB, Molecule* molecule);
-   void SetMessages();
-   double GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, 
-                             int mu, Molecule* molecule, double** gammaAB,
-                             double** orbitalElectronPopulation, double* atomicElectronPopulation,
-                             bool isGuess);
-   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);
-   void SetEnableAtomTypes();
 };
 
 ZindoS::ZindoS() : MolDS_cndo::Cndo2(){
@@ -55,6 +58,8 @@ void ZindoS::SetMessages(){
       = "Error in zindo::ZindoS::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in zindo::ZindoS::CheckEnableAtomType: Not enable 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->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,10 +79,9 @@ double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA
                                  double** orbitalElectronPopulation, double* atomicElectronPopulation,
                                  bool isGuess){
    double value;
-   value = atomA->GetIndoCoreIntegral(atomA->GetValence()[mu-firstAOIndexA], 
+   value = atomA->GetCoreIntegral(atomA->GetValence()[mu-firstAOIndexA], 
                                      gammaAB[atomAIndex][atomAIndex], 
-                                     isGuess);
-
+                                     isGuess, this->theory);
    if(!isGuess){
       double temp = 0.0;
       double coulomb = 0.0;
@@ -86,8 +90,8 @@ 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  = atomA->GetIndoCoulombInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex]);
-         exchange = atomA->GetIndoExchangeInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex]);
+         coulomb  = this->GetCoulombInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
+         exchange = this->GetExchangeInt(orbitalMu, orbitalLam, gammaAB[atomAIndex][atomAIndex], atomA);
          lammda = firstAOIndexA + v;
          temp += orbitalElectronPopulation[lammda][lammda]*(coulomb - 0.5*exchange);
       }
@@ -127,8 +131,8 @@ double ZindoS::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, i
       if(atomAIndex == atomBIndex){
          OrbitalType orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
          OrbitalType orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
-         coulomb  = atomA->GetIndoCoulombInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex]); 
-         exchange = atomA->GetIndoExchangeInt(orbitalMu, orbitalNu, gammaAB[atomAIndex][atomAIndex]); 
+         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];
       }
       else{
@@ -144,6 +148,76 @@ void ZindoS::CalcGammaAB(double** gammaAB, Molecule* molecule){
    // Do nothing;
 }
 
+// Apendix in [BZ_1972]
+// ZINDO Coulomb Interaction
+double ZindoS::GetCoulombInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+
+   double value=0.0;
+
+   // ToDo: Coulomb interaction
+   /*
+   if( orbital1 == s && orbital2 == s){ 
+      value = gamma;
+   }   
+   else if( orbital1 == s && ( orbital2 == px || orbital2 == py || orbital2 == pz )){ 
+      value = gamma;
+   }   
+   else if( (orbital1 == px || orbital1 == py || orbital1 == pz ) && orbital2 == s){ 
+      value = gamma;
+   }   
+   else if( (orbital1 == orbital2) && ( orbital1 == px || orbital1 == py || orbital1 == pz )){ 
+      value = gamma + 4.0*atom->GetIndoF2()/25.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;
+   }   
+   else{
+      cout << this->errorMessageCoulombInt;
+      cout << this->errorMessageAtomType << AtomTypeStr(atom->GetAtomType()) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
+      exit(EXIT_FAILURE);
+   }   
+   */
+   return value;
+
+}
+
+// Apendix in [BZ_1972]
+// ZINDO Exchange Interaction
+double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, double gamma, Atom* atom){
+
+   double value=0.0;
+
+   // ToDo: Exchange interaction
+   /*
+   if( orbital1 == orbital2){
+      value = this->GetCoulombInt(orbital1, orbital2, gamma, atom);
+   }   
+   else if( (orbital1 == s) && (orbital2 == px || orbital2 == py || orbital2 == pz ) ){
+      value = atom->GetIndoG1()/3.0;
+   }   
+   else if( (orbital1 == px || orbital1 == py || orbital1 == pz) && orbital2 == s  ){  
+      value = atom->GetIndoG1()/3.0;
+   }   
+   else if( (orbital1 != orbital2) 
+         && ( orbital1 == px || orbital1 == py || orbital1 == pz )
+         && ( orbital2 == px || orbital2 == py || orbital2 == pz ) ){
+      value = 3.0*atom->GetIndoF2()/25.0;
+   }   
+   else{
+      cout << this->errorMessageExchangeInt;
+      cout << this->errorMessageAtomType << AtomTypeStr(atom->GetAtomType()) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital1) << "\n";
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital2) << "\n";
+      exit(EXIT_FAILURE);
+   }   
+   */
+
+   return value;
+}
 
 }
 #endif