Commit: 5d7e889988c7cade162b71befb236cf1858b8598 - MolDS (git)

Commit MetaInfo

• Revision: 5d7e889988c7cade162b71befb236cf1858b8598 (tree)
• Zeit: 2010-12-08 22:32:05
• Autor: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Logic of core integrals for ZINDO is added.

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

Ändern Zusammenfassung

• modified: src/base/Parameters.h
• 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
• add: src/zindo/ZindoS.h

Diff

--- a/src/base/Parameters.h
+++ b/src/base/Parameters.h
@@ -23,6 +23,7 @@ private:
    void SetDefaultValues();
    double eV2AU;
    double angstrom2AU;
+   double kayser2AU;
    double bondingAdjustParameterK; //see (3.79) in J. A. Pople book
    TheoryType currentTheory;
 
@@ -36,6 +37,7 @@ public:
    void SetMaxIterationsSCF(int maxIterationsSCF);
    double GetEV2AU();
    double GetAngstrom2AU();
+   double GetKayser2AU();
    double GetBondingAdjustParameterK();
    TheoryType GetCurrentTheory();
    void SetCurrentTheory(TheoryType theory);
@@ -67,9 +69,9 @@ void Parameters::DeleteInstance(){
 void Parameters::SetDefaultValues(){
    this->eV2AU = 0.03674903;
    this->angstrom2AU = 1.0/0.5291772;
+   this->kayser2AU = 4.556336*pow(10.0,-6.0);
    this->thresholdSCF = pow(10.0, -8.0);
    this->maxIterationsSCF = 100;
-   this->bondingAdjustParameterK = 0.75;
    this->currentTheory = CNDO2;
 }
 
@@ -97,6 +99,10 @@ double Parameters::GetAngstrom2AU(){
    return this->angstrom2AU;
 }
 
+double Parameters::GetKayser2AU(){
+   return this->kayser2AU;
+}
+
 double Parameters::GetBondingAdjustParameterK(){
    return this->bondingAdjustParameterK;
 }

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -21,28 +21,52 @@ private:
    string errorMessageIndoExchangeInt;
    string errorMessageShellType;
    string errorMessageEffectivPrincipalQuantumNumber;
+   string errorMessageZindoCoreIntegral;
+   double GetJss();  // Part of Eq. (13) in [BZ_1979]
+   double GetJsp();  // Part of Eq. (13) in [BZ_1979]
+   double GetJsd();  // Part of Eq. (13) in [BZ_1979]
+   double GetJpp();  // Part of Eq. (13) in [BZ_1979]
+   double GetJpd();  // Part of Eq. (13) in [BZ_1979]
+   double GetJdd();  // Part of Eq. (13) in [BZ_1979]
 protected:
    string errorMessageIndoCoreIntegral;
+   string errorMessageZindoSCoreIntegral;
    string errorMessageAtomType;
    string errorMessageOrbitalType;
    double* xyz;
    AtomType atomType;
    vector<OrbitalType> valence;
-   double imuAmuS;
-   double imuAmuP;
-   double imuAmuD;
-   double bondingParameter;             // see Table 3.2 and 3.4 in J. A. Pople book
-   double coreCharge;                   // = Z_A
-   double effectiveNuclearChargeK;
-   double effectiveNuclearChargeL;
-   double effectiveNuclearChargeMsp;
-   double effectiveNuclearChargeMd;
    ShellType valenceShellType;
    int firstAOIndex;
-   int GetEffectivePrincipalQuantumNumber(ShellType shellType);
    int numberValenceElectrons;
-   double indoF2;                   // see (3.89) in J. A. Pople book
-   double indoG1;                   // see (3.88) in J. A. Pople book
+   double imuAmuS;                      // Table 3.4 or 3.5 in J. A. Pople book
+   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 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
+   double effectiveNuclearChargeMsp;    // Table 1.5 in J. A. Pople book
+   double effectiveNuclearChargeMd;     // Table 1.5 in J. A. Pople book
+   int GetEffectivePrincipalQuantumNumber(ShellType shellType); // Table 1.4 in J. A. Pople book
+   double indoF2;                   // (3.89) in J. A. Pople book
+   double indoG1;                   // (3.88) in J. A. Pople book
+   double zindoF0ss;                // Table 1 in ref. [RZ_1976], Table 1 in [AEZ_1986], or Table 1 in [GD_1972]
+   double zindoF0sd;                  // Table 1 in [AEZ_1986]
+   double zindoF0dd;                  // Table 1 in [AEZ_1986]
+   double zindoG1sp;                 // Table 3 in ref. [BZ_1979]
+   double zindoF2pp;                 // Table 3 in ref. [BZ_1979]
+   double zindoG2sd;                 // Table 3 in ref. [BZ_1979]
+   double zindoG1pd;                 // Table 3 in ref. [BZ_1979]
+   double zindoF2pd;                 // Table 3 in ref. [BZ_1979]
+   double zindoG3pd;                 // Table 3 in ref. [BZ_1979]
+   double zindoF2dd;                 // Table 3 in ref. [BZ_1979]
+   double zindoF4dd;                 // Table 3 in ref. [BZ_1979]
+   double IonPotS;   // Ionization potential, Table 4 in [BZ_1979]
+   double IonPotP;   // Ionization potential, Table 4 in [BZ_1979]
+   double IonPotD;   // Ionization potential, Table 4 in [BZ_1979]
+   double GetZindoCoreIntegral(OrbitalType orbital, int l, int m, int n); // Eq. (13) in [BZ_1979]
+
 public:
    Atom();
    Atom(double x, double y, double z);
@@ -61,7 +85,7 @@ public:
    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 GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess) = 0; // P82 - 83 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
 };
 
 Atom::Atom(){
@@ -89,12 +113,14 @@ void Atom::SetMessages(){
    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::GetIndoCoreIntegral: 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 = ";
    this->errorMessageOrbitalType = "\torbital type = ";
    this->errorMessageShellType = "\tshell type = ";
    this->errorMessageEffectivPrincipalQuantumNumber = 
       "Error in base::Atom::GetEffectivePrincipalQuantumNumber: invalid shelltype.\n";
+   this->errorMessageZindoCoreIntegral = "Error in base_stoms::Atom::GetZindoCoreINtegral: Invalid orbitalType.\n";
 }
 
 AtomType Atom::GetAtomType(){
@@ -276,6 +302,74 @@ double Atom::GetIndoExchangeInt(OrbitalType orbital1, OrbitalType orbital2, doub
    return value;
 }
 
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJss(){
+   return this->zindoF0ss;
+}
+
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJsp(){
+   // F0ss = F0sp
+   return this->zindoF0ss - this->zindoG1sp/6.0;
+}
+
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJsd(){
+   return this->zindoF0sd - this->zindoG2sd/10.0;
+}
+
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJpp(){
+   // F0pp = F0ss
+   return this->zindoF0ss - 2.0*this->zindoF2pp/25.0;
+}
+
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJpd(){
+   // F0pd = F0sd
+   return this->zindoF0sd - this->zindoG1pd/15.0 - 3.0*this->zindoG3pd/70.0;
+}
+
+// Part of Eq. (13) in [BZ_1979]
+double Atom::GetJdd(){
+   return this->zindoF0dd - 2.0*(this->zindoF2dd - this->zindoF4dd)/63.0;
+}
+
+// Eq. (13) in [BZ_1979]
+double Atom::GetZindoCoreIntegral(OrbitalType orbital, int l, int m, int n){
+   double value=0.0;
+
+   if(orbital == s){
+      value = -1.0*this->IonPotS 
+              - this->GetJss()*(double)(l-1) 
+              - this->GetJsp()*(double)m 
+              - this->GetJsd()*(double)n;
+   }
+   else if(orbital == px || orbital == py || orbital == pz){
+      value = -1.0*this->IonPotP
+              - this->GetJpp()*(double)(m-1) 
+              - this->GetJsp()*(double)l
+              - this->GetJpd()*(double)n;
+   }
+   else if(orbital == dxy || orbital == dyz || orbital == dzz || orbital == dzx || orbital == dxxyy ){
+      value = -1.0*this->IonPotD
+              - this->GetJdd()*(double)(n-1) 
+              - this->GetJsd()*(double)l
+              - this->GetJpd()*(double)m;
+   }
+   else{
+      cout << this->errorMessageZindoCoreIntegral;
+      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+      exit(EXIT_FAILURE);
+   }
+
+   return value;
+}
+
+
+
+
 }
 #endif
 

--- a/src/base/atoms/Catom.h
+++ b/src/base/atoms/Catom.h
@@ -11,7 +11,7 @@ class Catom : public Atom {
 private:
 public:
    Catom(double x, double y, double z);
-   double GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess); // P82 - 83 in J. A. Pople book.
+   double GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory); 
 };
 
 Catom::Catom(double x, double y, double z) : Atom(x, y, z){
@@ -33,31 +33,62 @@ Catom::Catom(double x, double y, double z) : Atom(x, y, z){
    this->numberValenceElectrons = 4;
    this->indoG1 = 0.267708;
    this->indoF2 = 0.17372;
+   this->zindoF0ss = 11.11 * Parameters::GetInstance()->GetEV2AU();                  
+   this->zindoF0sd = 0.0;                   
+   this->zindoF0dd = 0.0;                  
+   this->zindoG1sp = 55635*Parameters::GetInstance()->GetKayser2AU();                 
+   this->zindoF2pp = 36375*Parameters::GetInstance()->GetKayser2AU();                 
+   this->zindoG2sd = 0.0;                 
+   this->zindoG1pd = 0.0;                 
+   this->zindoF2pd = 0.0;                 
+   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->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }
 
-// P82 - 83 in J. A. Pople book.
-double Catom::GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess){
+double Catom::GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory){
    double value = 0.0;
-   if(orbital == s){
-      value = -1.0*this->imuAmuS;
-      if(!isGuess){
-         value -= (this->coreCharge-0.5)*gamma - (this->coreCharge - 1.5)*this->indoG1/6.0;
+
+   if(theory == INDO){
+      if(orbital == s){
+         value = -1.0*this->imuAmuS;
+         if(!isGuess){
+            value -= (this->coreCharge-0.5)*gamma - (this->coreCharge - 1.5)*this->indoG1/6.0;
+         }
       }
-   }
-   else if(orbital == px || orbital == py || orbital == pz){
-      value = -1.0*this->imuAmuP;
-      if(!isGuess){
-         value -= (this->coreCharge-0.5)*gamma 
-                 - this->indoG1/3.0 
-                 - (this->coreCharge - 2.5)*this->indoF2*2.0/25.0;
+      else if(orbital == px || orbital == py || orbital == pz){
+         value = -1.0*this->imuAmuP;
+         if(!isGuess){
+            value -= (this->coreCharge-0.5)*gamma 
+                    - this->indoG1/3.0 
+                    - (this->coreCharge - 2.5)*this->indoF2*2.0/25.0;
+         }
+      }
+      else{
+         cout << this->errorMessageIndoCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
       }
    }
-   else{
-      cout << this->errorMessageIndoCoreIntegral;
-      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
-      exit(EXIT_FAILURE);
+   else if(theory == ZINDOS){
+      if(orbital == s){
+         value = this->GetZindoCoreIntegral(orbital, 2, 2, 0);
+      }
+      else if(orbital == px || orbital == py || orbital == pz){
+         value = this->GetZindoCoreIntegral(orbital, 2, 2, 0);
+      }
+      else{
+         cout << this->errorMessageZindoSCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
+      }
    }
+
    return value;
 }
 

--- a/src/base/atoms/Hatom.h
+++ b/src/base/atoms/Hatom.h
@@ -12,7 +12,7 @@ class Hatom : public Atom {
 private:
 public:
    Hatom(double x, double y, double z);
-   double GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess); // P82 - 83 in J. A. Pople book.
+   double GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory); 
 };
 
 Hatom::Hatom(double x, double y, double z) : Atom(x, y, z){
@@ -31,23 +31,52 @@ Hatom::Hatom(double x, double y, double z) : Atom(x, y, z){
    this->numberValenceElectrons = 1;
    this->indoG1 = 0.0;
    this->indoF2 = 0.0;
+   this->zindoF0ss = 12.85 * Parameters::GetInstance()->GetEV2AU();                  
+   this->zindoF0sd = 0.0;                   
+   this->zindoF0dd = 0.0;                  
+   this->zindoG1sp = 0.0;                 
+   this->zindoF2pp = 0.0;                 
+   this->zindoG2sd = 0.0;                 
+   this->zindoG1pd = 0.0;                 
+   this->zindoF2pd = 0.0;                 
+   this->zindoG3pd = 0.0;                 
+   this->zindoF2dd = 0.0;                 
+   this->zindoF4dd = 0.0;                 
+   this->IonPotS = -13.06 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotP = 0.0 * Parameters::GetInstance()->GetEV2AU();
+   this->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }
 
-// P82 - 83 in J. A. Pople book.
-double Hatom::GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess){
+double Hatom::GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory){
    double value = 0.0;
-   if(orbital == s){
-      value = -1.0*this->imuAmuS;
-      if(!isGuess){
-         value -= 0.5*gamma;
+
+   if(theory == INDO){
+      if(orbital == s){
+         value = -1.0*this->imuAmuS;
+         if(!isGuess){
+            value -= 0.5*gamma;
+         }
+      }
+      else{
+         cout << this->errorMessageIndoCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
       }
    }
-   else{
-      cout << this->errorMessageIndoCoreIntegral;
-      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
-      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
-      exit(EXIT_FAILURE);
+   else if(theory == ZINDOS){
+      if(orbital == s){
+         value = this->GetZindoCoreIntegral(orbital, 1, 0, 0);
+      }
+      else{
+         cout << this->errorMessageZindoSCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
+      }
+      
    }
+
    return value;
 }
 

--- a/src/base/atoms/Liatom.h
+++ b/src/base/atoms/Liatom.h
@@ -12,7 +12,7 @@ class Liatom : public Atom {
 private:
 public:
    Liatom(double x, double y, double z);
-   double GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess); // P82 - 83 in J. A. Pople book.
+   double GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory); 
 };
 
 Liatom::Liatom(double x, double y, double z) : Atom(x, y, z){
@@ -34,30 +34,53 @@ Liatom::Liatom(double x, double y, double z) : Atom(x, y, z){
    this->numberValenceElectrons = 1;
    this->indoG1 = 0.092012;
    this->indoF2 = 0.049865;
+   this->zindoF0ss = 0.0;
+   this->zindoF0sd = 0.0;        
+   this->zindoF0dd = 0.0;              
+   this->zindoG1sp = 20194*Parameters::GetInstance()->GetKayser2AU();           
+   this->zindoF2pp = 10944*Parameters::GetInstance()->GetKayser2AU();        
+   this->zindoG2sd = 0.0;                
+   this->zindoG1pd = 0.0;                 
+   this->zindoF2pd = 0.0;                
+   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->IonPotD = 0.0 * Parameters::GetInstance()->GetEV2AU();
 }
 
-// P82 - 83 in J. A. Pople book.
-double Liatom::GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess){
+double Liatom::GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory){
    
    double value = 0.0;
-   if(orbital == s){
-      value = -1.0*this->imuAmuS;
-      if(!isGuess){
-         value -= 0.5*gamma;
+
+   if(theory == INDO){
+      if(orbital == s){
+         value = -1.0*this->imuAmuS;
+         if(!isGuess){
+            value -= 0.5*gamma;
+         }
       }
-   }
-   else if(orbital == px || orbital == py || orbital == pz){
-      value = -1.0*this->imuAmuP;
-      if(!isGuess){
-         value -= 0.5*gamma - this->indoG1/12.0;
+      else if(orbital == px || orbital == py || orbital == pz){
+         value = -1.0*this->imuAmuP;
+         if(!isGuess){
+            value -= 0.5*gamma - this->indoG1/12.0;
+         }
+      }
+      else{
+         cout << this->errorMessageIndoCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
       }
    }
-   else{
-      cout << this->errorMessageIndoCoreIntegral;
+   else if(theory == ZINDOS){
+      cout << this->errorMessageZindoSCoreIntegral;
       cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
       cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
       exit(EXIT_FAILURE);
    }
+
    return value;
 }
 

--- a/src/base/atoms/Satom.h
+++ b/src/base/atoms/Satom.h
@@ -11,7 +11,7 @@ class Satom : public Atom {
 private:
 public:
    Satom(double x, double y, double z);
-   double GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess); // P82 - 83 in J. A. Pople book.
+   double GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory); 
 };
 
 Satom::Satom(double x, double y, double z) : Atom(x, y, z){
@@ -39,17 +39,48 @@ Satom::Satom(double x, double y, double z) : Atom(x, y, z){
    this->effectiveNuclearChargeMsp = 5.45;
    this->effectiveNuclearChargeMd = 5.45;
    this->numberValenceElectrons = 6;
-   //this->indoG1 = 0.267708;
-   //this->indoF2 = 0.17372;
+   this->indoG1 = 0.267708;
+   this->indoF2 = 0.17372;
+   this->zindoF0ss = 8.96 * Parameters::GetInstance()->GetEV2AU();                  
+   this->zindoF0sd = 0.0;                   
+   this->zindoF0dd = 0.0;                 
+   this->zindoG1sp = 24807*Parameters::GetInstance()->GetKayser2AU();                 
+   this->zindoF2pp = 36600*Parameters::GetInstance()->GetKayser2AU();                 
+   this->zindoG2sd = 25972*Parameters::GetInstance()->GetKayser2AU();     
+   this->zindoG1pd = 34486*Parameters::GetInstance()->GetKayser2AU();        
+   this->zindoF2pd = 29173*Parameters::GetInstance()->GetKayser2AU();           
+   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();
 }
 
-// P82 - 83 in J. A. Pople book.
-double Satom::GetIndoCoreIntegral(OrbitalType orbital, double gamma, bool isGuess){
+double Satom::GetCoreIntegral(OrbitalType orbital, double gamma, bool isGuess, TheoryType theory){
    double value = 0.0;
-   cout << this->errorMessageIndoCoreIntegral;
-   cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
-   cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
-   exit(EXIT_FAILURE);
+
+   if(theory == INDO){
+      cout << this->errorMessageIndoCoreIntegral;
+      cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+      cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+      exit(EXIT_FAILURE);
+   }
+   else if(theory == ZINDOS){
+      if(orbital == s){
+         value = this->GetZindoCoreIntegral(orbital, 2, 4, 0);
+      }
+      else if(orbital == px || orbital == py || orbital == pz){
+         value = this->GetZindoCoreIntegral(orbital, 2, 4, 0);
+      }
+      else{
+         cout << this->errorMessageZindoSCoreIntegral;
+         cout << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         cout << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         exit(EXIT_FAILURE);
+      }
+   }
+
    return value;
 }
 

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -513,10 +513,12 @@ double Cndo2::GetFockOffDiagElement(Atom* atomA, Atom* atomB, int atomAIndex, in
                                     int mu, int nu, Molecule* molecule, double** gammaAB, double** overlap,
                                     double** orbitalElectronPopulation, bool isGuess){
    double value;
-   double K = 1.0;
+
+   double K = 1.0;  // = 1.0 or 0.75, see Eq. (3.79) in J. A. Pople book
    if(atomA->GetValenceShellType() >= m || atomB->GetValenceShellType() >= m){
       K = 0.75;
    }
+
    double bondParameter = 0.5*K*(atomA->GetBondingParameter() + atomB->GetBondingParameter()); 
    value =  bondParameter*overlap[mu][nu];
    if(!isGuess){

--- a/src/indo/Indo.h
+++ b/src/indo/Indo.h
@@ -76,9 +76,9 @@ double Indo::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;

--- /dev/null
+++ b/src/zindo/ZindoS.h
@@ -0,0 +1,152 @@
+#ifndef INCLUDED_INDO
+#define INCLUDED_INDO
+
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include<iostream>
+#include<vector>
+#include"../cndo/Cndo2.h"
+
+using namespace std;
+using namespace MolDS_base;
+using namespace MolDS_base_atoms;
+
+namespace MolDS_zindo{
+
+/***
+ *  Refferences for Indo are [PB_1970] and [PS_1966].
+ */
+class ZindoS : public MolDS_cndo::Cndo2{
+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(){
+   this->theory = ZINDOS;
+   this->SetMessages();
+   this->SetEnableAtomTypes();
+   //cout << "ZindoS created\n";
+}
+
+ZindoS::~ZindoS(){
+   //cout << "ZindoS deleted\n";
+}
+
+void ZindoS::SetMessages(){
+   this->errorMessageSCFNotConverged 
+      = "Error in zindo::ZindoS::DoesSCF: SCF did not met convergence criterion. maxIterationsSCF=";
+   this->errorMessageMoleculeNotSet 
+      = "Error in zindo::ZindoS::DoesSCF: A molecule is not set.\n";
+   this->errorMessageOddTotalValenceElectrions 
+      = "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->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";
+}
+
+void ZindoS::SetEnableAtomTypes(){
+   this->enableAtomTypes.clear();
+   this->enableAtomTypes.push_back(H);
+   this->enableAtomTypes.push_back(C);
+   this->enableAtomTypes.push_back(N);
+   this->enableAtomTypes.push_back(O);
+   this->enableAtomTypes.push_back(S);
+}
+
+double ZindoS::GetFockDiagElement(Atom* atomA, int atomAIndex, int firstAOIndexA, int mu, 
+                                 Molecule* molecule, double** gammaAB,
+                                 double** orbitalElectronPopulation, double* atomicElectronPopulation,
+                                 bool isGuess){
+   double value;
+   value = atomA->GetIndoCoreIntegral(atomA->GetValence()[mu-firstAOIndexA], 
+                                     gammaAB[atomAIndex][atomAIndex], 
+                                     isGuess);
+
+   if(!isGuess){
+      double temp = 0.0;
+      double coulomb = 0.0;
+      double exchange = 0.0;
+      int lammda = 0;
+      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]);
+         lammda = firstAOIndexA + v;
+         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];
+         }
+      }
+      value += temp;
+   }
+
+   return value;
+}
+
+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()); 
+
+   if(isGuess){
+      value = bondParameter*overlap[mu][nu];
+   }
+   else{
+      double coulomb = 0.0;
+      double exchange = 0.0;
+      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]); 
+         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];
+      }
+   }
+
+   return value;
+}
+
+void ZindoS::CalcGammaAB(double** gammaAB, Molecule* molecule){
+   // Do nothing;
+}
+
+
+}
+#endif
+
+
+