コミット: a8a0f2ebad34e652113c80007c98e0f86dabe55a - MolDS (git)

コミットメタ情報

• リビジョン: a8a0f2ebad34e652113c80007c98e0f86dabe55a (tree)
• 日時: 2012-11-07 19:37:22
• 作者: Katsuhiko Nishimra <ktns.87@gmai...>
• コミッター: Katsuhiko Nishimra

ログメッセージ

Merge trunk into branch gdiis. #28915

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/branches/gdiis@1104 1136aad2-a195-0410-b898-f5ea1d11b9d8

変更サマリ

• modified: src/base/ElectronicStructure.h
• modified: src/base/InputParser.cpp
• modified: src/base/InputParser.h
• modified: src/base/MathUtilities.cpp
• modified: src/base/MathUtilities.h
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/indo/Indo.cpp
• modified: src/nasco/NASCO.cpp
• modified: src/nasco/NASCO.h
• modified: src/optimization/BFGS.cpp
• modified: src/wrappers/Blas.cpp
• modified: src/wrappers/Blas.h
• modified: src/wrappers/Blas_GNU.cpp
• modified: src/wrappers/Lapack.cpp
• modified: src/wrappers/Lapack.h
• modified: src/wrappers/Lapack_GNU.cpp
• modified: src/zindo/ZindoS.cpp
• modified: src/zindo/ZindoS.h

差分

--- a/src/base/ElectronicStructure.h
+++ b/src/base/ElectronicStructure.h
@@ -43,6 +43,12 @@ public:
    virtual void CalcOverlapMOsWithAnotherElectronicStructure(double** overlapMOs, 
                                                              double const* const* overlapAOs,
                                                              const MolDS_base::ElectronicStructure& lhsElectronicStructure) const = 0;
+   virtual void CalcOverlapSingletSDsWithAnotherElectronicStructure(double** overlapSingletSDs, 
+                                                                    double const* const* overlapMOs) const = 0;
+   virtual void CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs, 
+                                                             double const* const* overlapSingletSDs,
+                                                             const MolDS_base::ElectronicStructure& lhsElectronicStructure) const = 0;
+
 };
 
 }

--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2011-2012 Michihiro Okuyama                              // 
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -72,6 +73,8 @@ void InputParser::DeleteInstance(){
 }
 
 void InputParser::SetMessages(){
+   this->errorMessageInputFileEmpty
+      = "Error in base::InputParser::GetInputTerms: Input file is empty.\n"; 
    this->errorMessageNotFoundInputFile
       = "Error in base::InputParser::StoreInputTermsFromFile: Not found.\n"; 
    this->errorMessageNonValidExcitedStatesMD
@@ -368,6 +371,15 @@ vector<string> InputParser::GetInputTerms(int argc, char *argv[]) const{
       char* fileName = argv[1];
       this->StoreInputTermsFromFile(inputTerms,fileName);
    }
+   if (inputTerms.size() == 0) {
+      char* fileName = argv[1];
+      stringstream ss; 
+      ss << this->errorMessageInputFileEmpty;  
+      if (argc > 1) {
+         ss << this->errorMessageInputFile << fileName << endl; 
+      }
+      throw MolDSException(ss.str()); 
+   }
    return inputTerms;
 }
 

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2011-2012 Michihiro Okuyama                              // 
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -31,6 +32,7 @@ private:
    InputParser();
    ~InputParser();
    void SetMessages();
+   std::string errorMessageInputFileEmpty;
    std::string errorMessageNotFoundInputFile; 
    std::string errorMessageNonValidExcitedStatesMD;
    std::string errorMessageNonValidExcitedStatesMC;

--- a/src/base/MathUtilities.cpp
+++ b/src/base/MathUtilities.cpp
@@ -22,9 +22,14 @@
 #include<sstream>
 #include<math.h>
 #include<stdexcept>
+#include<boost/format.hpp>
+#include"PrintController.h"
 #include"MolDSException.h"
+#include"Uncopyable.h"
+#include"../wrappers/Lapack.h"
 #include"Enums.h"
 #include"MathUtilities.h"
+#include"MallocerFreer.h"
 using namespace std;
 
 namespace MolDS_base{
@@ -133,5 +138,21 @@ void CalcRotatingMatrix(double matrix[][3], double sita, CartesianType cartesian
       throw MolDSException(ss.str());
    }
 }
+
+// calculate determinant of the matrix
+double GetDeterminant(double** matrix, int dim){
+   double determinant=1.0;
+   int* ipiv=NULL;
+   MallocerFreer::GetInstance()->Malloc<int>(&ipiv, dim);
+   MolDS_wrappers::Lapack::GetInstance()->Dgetrf(matrix, ipiv, dim, dim);
+   for(int i=0; i<dim; i++){
+      determinant*=matrix[i][i];
+      if(ipiv[i] != i-1){
+         determinant *= -1.0;
+      }
+   }
+   MallocerFreer::GetInstance()->Free<int>(&ipiv, dim);
+   return determinant;
+}
 }
  

--- a/src/base/MathUtilities.h
+++ b/src/base/MathUtilities.h
@@ -29,6 +29,8 @@ template <typename T> T Max(T a, T b);
 template <typename T> T min(T a, T b);
 // rotating matrix
 void CalcRotatingMatrix(double matrix[][3], double sita, CartesianType cartesianType);
+// calculate determinant of the matrix. Note taht the matrix will be destroid
+double GetDeterminant(double** matrix, int dim);
 }
 #endif
  

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -184,6 +184,8 @@ void Cndo2::SetMessages(){
       = "Error in cndo::Cndo2::CalcOverlapAOsDifferentConfigurations: Number Atoms in lhs and rhs are different.\n";
    this->errorMessageCalcOverlapAOsDifferentConfigurationsOverlapAOsNULL
       = "Error in cndo::Cndo2::CalcOverlapAOsDifferentConfigurations: ovelrapAOs is NULL.\n";
+   this->errorMessageNonExcitedStates 
+      = "Error in cndo::CNDO2::Excited states can not be calculated with CNDO2.\n";
    this->errorMessageLhs = "lhs: ";
    this->errorMessageRhs = "rhs: ";
    this->errorMessageFromState = "\tfrom state = ";
@@ -3580,6 +3582,8 @@ void Cndo2::CalcOverlapMOsWithAnotherElectronicStructure(double** overlapMOs,
    double const* const* rhsFockMatrix = this->fockMatrix;
    double const* const* lhsFockMatrix = lhsElectronicStructure.GetFockMatrix();
    int totalAONumber = this->molecule->GetTotalNumberAOs();
+   int usedMONumber = this->molecule->GetTotalNumberValenceElectrons()/2
+                     +Parameters::GetInstance()->GetActiveVirCIS();
    double** tmpMatrix=NULL;
    try{
       MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,totalAONumber,totalAONumber);
@@ -3589,13 +3593,13 @@ void Cndo2::CalcOverlapMOsWithAnotherElectronicStructure(double** overlapMOs,
       double beta=0.0;
       MolDS_wrappers::Blas::GetInstance()->Dgemm(isColumnMajorOverlapAOs,
                                                  isColumnMajorRhsFock,
-                                                 totalAONumber,totalAONumber,totalAONumber,
+                                                 totalAONumber,usedMONumber,totalAONumber,
                                                  alpha,
                                                  overlapAOs,
                                                  rhsFockMatrix,
                                                  beta,
                                                  tmpMatrix);
-      MolDS_wrappers::Blas::GetInstance()->Dgemm(totalAONumber,totalAONumber,totalAONumber,
+      MolDS_wrappers::Blas::GetInstance()->Dgemm(usedMONumber,totalAONumber,totalAONumber,
                                                  lhsFockMatrix,
                                                  tmpMatrix,
                                                  overlapMOs);
@@ -3608,6 +3612,33 @@ void Cndo2::CalcOverlapMOsWithAnotherElectronicStructure(double** overlapMOs,
    MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,totalAONumber,totalAONumber);
 }
 
+// calculate OverlapSingletSDs matrix between different electronic-structure, S^{SSD}_{ij}.
+// i and j are singlet SDs belonging to left and right hand side electronic-structures, respectively.
+// The index i=0 means the Hartree-Fock state.
+// This overlapsingletSDs are calculated from overlapMOs.
+// Note that rhs-electronic-structure is this electronic-structure  
+// and lhs-electronic-structure is another electronic-structure.
+void Cndo2::CalcOverlapSingletSDsWithAnotherElectronicStructure(double** overlapSingletSDs, 
+                                                                double const* const* overlapMOs) const{
+      stringstream ss;
+      ss << this->errorMessageNonExcitedStates;
+      throw MolDSException(ss.str());
+}
+
+// calculate overlapESs (ES means eigenstate) matrix between different electronic-structure, S^{ES}_{ij}.
+// i and j are singlet SDs belonging to left and right hand side electronic-structures, respectively.
+// The index i=0 means the ground state.
+// This overlapESs is calculated from the overlapsingletSDs.
+// Note that rhs-electronic-structure is this electronic-structure  
+// and lhs-electronic-structure is another electronic-structure.
+void Cndo2::CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs, 
+                                                         double const* const* overlapSingletSDs,
+                                                         const MolDS_base::ElectronicStructure& lhsElectronicStructure) const{
+      stringstream ss;
+      ss << this->errorMessageNonExcitedStates;
+      throw MolDSException(ss.str());
+}
+
 // calculate OverlapAOs matrix. E.g. S_{\mu\nu} in (3.74) in J. A. Pople book.
 void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
    int totalAONumber = molecule.GetTotalNumberAOs();

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -46,6 +46,11 @@ public:
    void CalcOverlapMOsWithAnotherElectronicStructure(double** overlapMOs,
                                                      double const* const* overlapAOs,
                                                      const MolDS_base::ElectronicStructure& lhsElectronicStructure) const;
+   virtual void CalcOverlapSingletSDsWithAnotherElectronicStructure(double** overlapSingletSDs, 
+                                                                    double const* const* overlapMOs) const;
+   virtual void CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs, 
+                                                             double const* const* overlapSingletSDs,
+                                                             const MolDS_base::ElectronicStructure& lhsElectronicStructure) const;
    MolDS_base::TheoryType GetTheoryType() const;
 protected:
    std::string errorMessageAtomA;
@@ -74,6 +79,7 @@ protected:
    std::string errorMessageCalcFrequenciesNormalModesBadTheory;
    std::string errorMessageFromState;
    std::string errorMessageToState;
+   std::string errorMessageNonExcitedStates;
    std::string messageSCFMetConvergence;
    std::string messageStartSCF;
    std::string messageDoneSCF;

--- a/src/indo/Indo.cpp
+++ b/src/indo/Indo.cpp
@@ -82,6 +82,8 @@ void Indo::SetMessages(){
       = "Error in indo::Indo::GetElectronicEnergy: excitedEnergies is NULL\n";
    this->errorMessageCalcFrequenciesNormalModesBadTheory
       = "Error in indo::Indo::CalcFrequenciesNormalModesBadTheory: INDO is not supported for frequency (normal mode) analysis.\n";
+   this->errorMessageNonExcitedStates 
+      = "Error in indo::Indo::Excited states can not be calculated with INDO.\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";

--- a/src/nasco/NASCO.cpp
+++ b/src/nasco/NASCO.cpp
@@ -105,13 +105,14 @@ void NASCO::DoNASCO(Molecule& molecule){
    molecule.OutputXyzCOC();
    molecule.OutputMomenta();
 
-   // malloc ovelap AOs, MOs, ESs between differentMolecules
+   // malloc ovelap AOs, MOs, Singlet Slater Determinants, and Eigenstates between differentMolecules
    double** overlapAOs = NULL;
    double** overlapMOs = NULL;
+   double** overlapSingletSDs = NULL;
    double** overlapESs = NULL;
 
    try{
-      this->MallocOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapESs, molecule);
+      this->MallocOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapSingletSDs, &overlapESs, molecule);
       for(int s=0; s<totalSteps; s++){
          this->OutputLog(boost::format("%s%d\n") % this->messageStartStepNASCO.c_str() % (s+1) );
 
@@ -141,7 +142,7 @@ void NASCO::DoNASCO(Molecule& molecule){
          // update momenta
          this->UpdateMomenta(molecule, matrixForce, dt);
 
-         // calculate overlaps
+         // calculate overlaps 
          currentES->CalcOverlapAOsWithAnotherConfiguration(overlapAOs, tmpMolecule);
          cout << "overlapAOs" << endl;
          for(int i=0; i<molecule.GetTotalNumberAOs(); i++){
@@ -160,6 +161,31 @@ void NASCO::DoNASCO(Molecule& molecule){
             }
             cout << endl;
          }
+         cout << endl;
+
+         currentES->CalcOverlapSingletSDsWithAnotherElectronicStructure(overlapSingletSDs, overlapMOs);
+         int dimOverlapSingletSDs = Parameters::GetInstance()->GetActiveOccCIS()
+                                  *Parameters::GetInstance()->GetActiveOccCIS()
+                                  +1;
+         cout << "overlap singlet slater determinants" << endl;
+         for(int i=0; i<dimOverlapSingletSDs; i++){
+            for(int j=0; j<dimOverlapSingletSDs; j++){
+               printf("%e\t",overlapSingletSDs[i][j]);
+            }
+            cout << endl;
+         }
+         cout << endl;
+
+         currentES->CalcOverlapESsWithAnotherElectronicStructure(overlapESs, overlapSingletSDs, *tmpES);
+         int dimOverlapESs = Parameters::GetInstance()->GetNumberElectronicStatesNASCO();
+         cout << "overlapESs" << endl;
+         for(int i=0; i<dimOverlapESs; i++){
+            for(int j=0; j<dimOverlapESs; j++){
+               printf("%e\t",overlapESs[i][j]);
+            }
+            cout << endl;
+         }
+         cout << endl;
             
          // Synchronous molecular configuration and electronic states
          this->SynchronousMolecularConfiguration(molecule, tmpMolecule);
@@ -179,10 +205,10 @@ void NASCO::DoNASCO(Molecule& molecule){
       }
    }
    catch(MolDSException ex){
-      this->FreeOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapESs, molecule);
+      this->FreeOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapSingletSDs, &overlapESs, molecule);
       throw ex;
    }
-   this->FreeOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapESs, molecule);
+   this->FreeOverlapsDifferentMolecules(&overlapAOs, &overlapMOs, &overlapSingletSDs, &overlapESs, molecule);
    this->OutputLog(this->messageEndNASCO);
 }
 
@@ -305,26 +331,36 @@ void NASCO::CheckEnableTheoryType(TheoryType theoryType){
 
 void NASCO::MallocOverlapsDifferentMolecules(double*** overlapAOs,
                                              double*** overlapMOs, 
+                                             double*** overlapSingletSDs, 
                                              double*** overlapESs, 
                                              const Molecule& molecule) const{
    int dimOverlapAOs = molecule.GetTotalNumberAOs();
    int dimOverlapMOs = dimOverlapAOs;
+   int dimOverlapSingletSDs = Parameters::GetInstance()->GetActiveOccCIS()
+                            *Parameters::GetInstance()->GetActiveVirCIS()
+                            +1;
    int dimOverlapESs = Parameters::GetInstance()->GetNumberElectronicStatesNASCO();
-   MallocerFreer::GetInstance()->Malloc<double>(overlapAOs, dimOverlapAOs, dimOverlapAOs);
-   MallocerFreer::GetInstance()->Malloc<double>(overlapMOs, dimOverlapMOs, dimOverlapMOs);
-   MallocerFreer::GetInstance()->Malloc<double>(overlapESs, dimOverlapESs, dimOverlapESs);
+   MallocerFreer::GetInstance()->Malloc<double>(overlapAOs,        dimOverlapAOs,        dimOverlapAOs);
+   MallocerFreer::GetInstance()->Malloc<double>(overlapMOs,        dimOverlapMOs,        dimOverlapMOs);
+   MallocerFreer::GetInstance()->Malloc<double>(overlapSingletSDs, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Malloc<double>(overlapESs,        dimOverlapESs,        dimOverlapESs);
 }
 
 void NASCO::FreeOverlapsDifferentMolecules(double*** overlapAOs,
                                            double*** overlapMOs, 
+                                           double*** overlapSingletSDs, 
                                            double*** overlapESs, 
                                            const MolDS_base::Molecule& molecule) const{
    int dimOverlapAOs = molecule.GetTotalNumberAOs();
    int dimOverlapMOs = dimOverlapAOs;
+   int dimOverlapSingletSDs = Parameters::GetInstance()->GetActiveOccCIS()
+                            *Parameters::GetInstance()->GetActiveVirCIS()
+                            +1;
    int dimOverlapESs = Parameters::GetInstance()->GetNumberElectronicStatesNASCO();
-   MallocerFreer::GetInstance()->Free<double>(overlapAOs, dimOverlapAOs, dimOverlapAOs);
-   MallocerFreer::GetInstance()->Free<double>(overlapMOs, dimOverlapMOs, dimOverlapMOs);
-   MallocerFreer::GetInstance()->Free<double>(overlapESs, dimOverlapESs, dimOverlapESs);
+   MallocerFreer::GetInstance()->Free<double>(overlapAOs,        dimOverlapAOs,        dimOverlapAOs);
+   MallocerFreer::GetInstance()->Free<double>(overlapMOs,        dimOverlapMOs,        dimOverlapMOs);
+   MallocerFreer::GetInstance()->Free<double>(overlapSingletSDs, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Free<double>(overlapESs,        dimOverlapESs,        dimOverlapESs);
 }
 
 }

--- a/src/nasco/NASCO.h
+++ b/src/nasco/NASCO.h
@@ -58,10 +58,12 @@ private:
    double OutputEnergies(const MolDS_base::ElectronicStructure& electronicStructure, const MolDS_base::Molecule& molecule);
    void MallocOverlapsDifferentMolecules(double*** overlapAOs,
                                          double*** overlapMOs, 
+                                         double*** overlapSingleSDs,
                                          double*** overlapESs, 
                                          const MolDS_base::Molecule& molecule) const;
    void FreeOverlapsDifferentMolecules(double*** overlapAOs,
                                        double*** overlapMOs, 
+                                       double*** overlapSingleSDs,
                                        double*** overlapESs, 
                                        const MolDS_base::Molecule& molecule) const;
 };

--- a/src/optimization/BFGS.cpp
+++ b/src/optimization/BFGS.cpp
@@ -32,6 +32,7 @@
 #include"../base/PrintController.h"
 #include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
+#include"../wrappers/Blas.h"
 #include"../wrappers/Lapack.h"
 #include"../base/Enums.h"
 #include"../base/MallocerFreer.h"
@@ -357,90 +358,40 @@ void BFGS::UpdateHessian(double **matrixHessian,
                          double const* vectorOldForce,
                          double const* vectorDisplacement) const{
    double const* const K = &vectorDisplacement[0]; // K_k in eq. 15 on [SJTO_1983]
-   double*  P    = NULL;                           // P_k in eq. 14 on [SJTO_1983]
-   double** PP   = NULL;                           // P_k P_k^T at second term in RHS of Eq. (13) in [SJTO_1983]
-   double*  HK   = NULL;                           // H_k K_k at third term on RHS of Eq. (13) in [SJTO_1983]
-   double** HKKH = NULL;                           // H_k K_k K_k^T H_k at third term on RHS of Eq. (13) in [SJTO_1983]
+   double *P  = NULL;                              // P_k in eq. 14 on [SJTO_1983]
+   double *HK = NULL;                              // H_k K_k at third term on RHS of Eq. (13) in [SJTO_1983]
    try{
       MallocerFreer::GetInstance()->Malloc(&P, dimension);
-      MallocerFreer::GetInstance()->Malloc(&PP, dimension, dimension);
       MallocerFreer::GetInstance()->Malloc(&HK, dimension);
-      MallocerFreer::GetInstance()->Malloc(&HKKH, dimension,dimension);
       double KHK = 0;
       double PK = 0;
-#pragma omp parallel for schedule(auto)
-      for(int i=0; i<dimension; i++){
-         // initialize P_k according to Eq. (14) in [SJTO_1983]
-         // note: gradient = -1 * force
-         P[i] = - (vectorForce[i] - vectorOldForce[i]);
-      }
+      // initialize P_k according to Eq. (14) in [SJTO_1983]
+      // note: gradient = -1 * force
+      MolDS_wrappers::Blas::GetInstance()->Dcopy(dimension, vectorOldForce, P);
+      MolDS_wrappers::Blas::GetInstance()->Daxpy(dimension, -1.0, vectorForce, P);
 
       // P_k^T K_k at second term at RHS of Eq. (13) in [SJTO_1983]
-#pragma omp parallel
-      {
-#pragma omp for schedule(auto) reduction(+:PK)
-         for(int i=0; i<dimension;i++){
-            PK += P[i] * K[i];
-         }
-
-         //P_k P_k^T at second term in RHS of Eq. (13) in [SJTO_1983]
-         for(int i=0; i<dimension;i++){
-#pragma omp for schedule(auto)
-            for(int j=i;j<dimension;j++){
-               PP[i][j] = PP[j][i] = P[i] * P[j];
-            }
-         }
+      PK = MolDS_wrappers::Blas::GetInstance()->Ddot(dimension, P, K);
 
-         //H_k K_k at third term on RHS of Eq. (13) in [SJTO_1983]
-#pragma omp for schedule(auto)
-         for(int i=0; i<dimension; i++){
-            double tmp=0;
-            for(int j=0; j<dimension; j++){
-               tmp += matrixHessian[i][j] * K[j];
-            }
-            HK[i] = tmp;
-         }
-      }
+      //H_k K_k at third term on RHS of Eq. (13) in [SJTO_1983]
+      MolDS_wrappers::Blas::GetInstance()->Dsymv(dimension, matrixHessian, K, HK);
 
       //K_k^T H_k K_k at third term on RHS of Eq. (13) in [SJTO_1983]
-#pragma omp parallel
-      {
-#pragma omp for schedule(auto) reduction(+:KHK)
-         for(int i=0;i<dimension;i++){
-            KHK += K[i]*HK[i];
-         }
-
-         //H_k K_k K_k^T H_k at third term on RHS of Eq. (13) in [SJTO_1983]
-         for(int i=0;i<dimension;i++){
-#pragma omp for schedule(auto)
-            for(int j=i;j<dimension;j++){
-               // H_k K_k = (K_k^T H_k)^T because H_k^T = H_k
-               HKKH[i][j] = HKKH[j][i] = HK[i] * HK[j];
-            }
-         }
-      }
+      KHK = MolDS_wrappers::Blas::GetInstance()->Ddot(dimension, K, HK);
 
       // Calculate H_k+1 according to Eq. (13) in [SJTO_1983]
-      for(int i=0;i<dimension;i++){
-#pragma omp parallel for schedule(auto)
-         for(int j=i;j<dimension;j++){
-            matrixHessian[i][j]+= (PP[i][j] / PK
-                                - HKKH[i][j] / KHK);
-            matrixHessian[j][i] = matrixHessian[i][j];
-         }
-      }
+      // Add second term in RHS of Eq. (13) in [SJTO_1983]
+      MolDS_wrappers::Blas::GetInstance()->Dsyr(dimension,  1.0/PK,  P,  matrixHessian);
+      // Add third term in RHS of Eq. (13) in [SJTO_1983]
+      MolDS_wrappers::Blas::GetInstance()->Dsyr(dimension, -1.0/KHK, HK, matrixHessian);
    }
    catch(MolDSException ex){
       MallocerFreer::GetInstance()->Free(&P, dimension);
-      MallocerFreer::GetInstance()->Free(&PP, dimension, dimension);
       MallocerFreer::GetInstance()->Free(&HK, dimension);
-      MallocerFreer::GetInstance()->Free(&HKKH, dimension,dimension);
       throw ex;
    }
    MallocerFreer::GetInstance()->Free(&P, dimension);
-   MallocerFreer::GetInstance()->Free(&PP, dimension, dimension);
    MallocerFreer::GetInstance()->Free(&HK, dimension);
-   MallocerFreer::GetInstance()->Free(&HKKH, dimension,dimension);
 }
 
 // Level shift eigenvalues of redandant modes to largeEigenvalue

--- a/src/wrappers/Blas.cpp
+++ b/src/wrappers/Blas.cpp
@@ -57,6 +57,66 @@ void Blas::DeleteInstance(){
    blas = NULL;
 }
 
+// vectorY = vectorX
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dcopy(int n,
+                 double const* vectorX,
+                 double *      vectorY)const{
+   int incrementX=1;
+   int incrementY=1;
+   this->Dcopy(n, vectorX, incrementX, vectorY, incrementY);
+}
+
+// vectorY = vectorX
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dcopy(int n,
+                 double const* vectorX, int incrementX,
+                 double*       vectorY, int incrementY) const{
+   dcopy(&n, vectorX, &incrementX, vectorY, &incrementY);
+}
+
+// vectorY = alpha*vectorX + vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Daxpy(int n, double alpha,
+           double const* vectorX,
+           double*       vectorY) const{
+   int incrementX=1;
+   int incrementY=1;
+   this->Daxpy(n, alpha, vectorX, incrementX, vectorY, incrementY);
+}
+
+// vectorY = alpha*vectorX + vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Daxpy(int n, double alpha,
+           double const* vectorX, int incrementX,
+           double*       vectorY, int incrementY) const{
+   daxpy(&n, &alpha, vectorX, &incrementX, vectorY, &incrementY);
+}
+
+// returns vectorX^T*vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+double Blas::Ddot(int n,
+            double const* vectorX,
+            double const* vectorY) const{
+   int incrementX=1;
+   int incrementY=1;
+   return this->Ddot(n, vectorX, incrementX, vectorY, incrementY);
+}
+
+// returns vectorX^T*vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+double Blas::Ddot(int n,
+            double const* vectorX, int incrementX,
+            double const* vectorY, int incrementY)const{
+   return ddot(&n, vectorX, &incrementX, vectorY, &incrementY);
+}
+
 // vectorY = matrixA*vectorX
 //    matrixA: m*n-matrix (matrixA[m][n] in row-major (C/C++ style))
 //    vectorX: n-vector
@@ -74,7 +134,7 @@ void Blas::Dgemv(int m, int n,
 }
 
 // vectorY = alpha*matrixA*vectorX + beta*vectorY
-//    matrixA: m*n-matrix (matrixA[m][n] in row-major (C/C++ style))
+//    matrixA: m*n-matrix
 //    vectorX: n-vector
 //    vectorY: m-vector
 void Blas::Dgemv(bool isColumnMajorMatrixA,
@@ -99,6 +159,60 @@ void Blas::Dgemv(bool isColumnMajorMatrixA,
    dgemv(&transA, &m, &n, &alpha, a, &lda, vectorX, &incrementX, &beta, vectorY, &incrementY);
 }
 
+// vectorY = matrixA*vectorX
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part)
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dsymv(int n,
+           double const* const* matrixA,
+           double const* vectorX,
+           double*       vectorY) const{
+   bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major.
+   int incrementX=1, incrementY=1;
+   double alpha=1.0, beta=0.0;
+   this->Dsymv(n, alpha, matrixA, vectorX, incrementX, beta, vectorY, incrementY);
+}
+
+// vectorY = alpha*matrixA*vectorX + beta*vectorY
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part)
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dsymv(int n, double alpha,
+           double const* const* matrixA,
+           double const* vectorX, int incrementX,
+           double beta,
+           double*       vectorY, int incrementY) const{
+   double const* a = &matrixA[0][0];
+   char uploA='U';
+   int lda = n;
+   dsymv(&uploA, &n, &alpha, a, &lda, vectorX, &incrementX, &beta, vectorY, &incrementY);
+}
+
+// matrixA = alpha*vectorX*vectorX^T + matrixA
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.)
+//    vectorX: n-matrix
+void Blas::Dsyr(int n, double alpha,
+          double const* vectorX,
+          double ** matrixA)const{
+   int incrementX=1;
+   this->Dsyr(n, alpha, vectorX, incrementX, matrixA);
+}
+
+void Blas::Dsyr(int n, double alpha,
+          double const* vectorX, int incrementX,
+          double ** matrixA)const{
+   double* a = &matrixA[0][0];
+   char uploA='U';
+   int lda = n;
+   dsyr(&uploA, &n, &alpha, vectorX, &incrementX, a, &lda);
+#pragma omp parallel for schedule(auto)
+   for(int i=0;i<n;i++){
+      for(int j=i+1;j<n;j++){
+         matrixA[j][i] = matrixA[i][j];
+      }
+   }
+}
+
 // matrixC = matrixA*matrixB
 //    matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style))
 //    matrixB: k*n-matrix (matrixB[k][n] in row-major (C/C++ style))
@@ -115,8 +229,8 @@ void Blas::Dgemm(int m, int n, int k,
 }
 
 // matrixC = alpha*matrixA*matrixB + beta*matrixC
-//    matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style))
-//    matrixB: k*n-matrix (matrixB[k][n] in row-major (C/C++ style))
+//    matrixA: m*k-matrix 
+//    matrixB: k*n-matrix
 //    matrixC: m*n-matrix (matrixC[m][n] in row-major (C/C++ style))
 void Blas::Dgemm(bool isColumnMajorMatrixA, 
                  bool isColumnMajorMatrixB, 

--- a/src/wrappers/Blas.h
+++ b/src/wrappers/Blas.h
@@ -24,6 +24,24 @@ class Blas: public MolDS_base::PrintController, private MolDS_base::Uncopyable{
 public:
    static Blas* GetInstance();
    static void DeleteInstance();
+   void Dcopy(int n,
+              double const* vectorX,
+              double*       vectorY) const;
+   void Dcopy(int n,
+              double const* vectorX, int incrementX,
+              double*       vectorY, int incrementY) const;
+   void Daxpy(int n, double alpha,
+              double const* vectorX,
+              double*       vectorY) const;
+   void Daxpy(int n, double alpha,
+              double const* vectorX, int incrementX,
+              double*       vectorY, int incrementY) const;
+   double Ddot(int n,
+               double const* vectorX,
+               double const* vectorY) const;
+   double Ddot(int n,
+               double const* vectorX, int incrementX,
+               double const* vectorY, int incrementY)const;
    void Dgemv(int m, int n,
               double const* const* matrixA,
               double const* vectorX,
@@ -37,6 +55,21 @@ public:
               double beta,
               double* vectorY,
               int incrementY) const;
+   void Dsymv(int n,
+              double const* const* matrixA,
+              double const* vectorX,
+              double* vectorY) const;
+   void Dsymv(int n, double alpha,
+              double const* const* matrixA,
+              double const* vectorX, int incrementX,
+              double beta,
+              double*       vectorY, int incrementY) const;
+   void Dsyr(int n, double alpha,
+             double const* vectorX,
+             double ** matrixA)const;
+   void Dsyr(int n, double alpha,
+             double const* vectorX, int incrementX,
+             double ** matrixA)const;
    void Dgemm(int m, int n, int k, 
               double const* const* matrixA, 
               double const* const* matrixB, 

--- a/src/wrappers/Blas_GNU.cpp
+++ b/src/wrappers/Blas_GNU.cpp
@@ -57,6 +57,70 @@ void Blas::DeleteInstance(){
    blas = NULL;
 }
 
+// vectorY = vectorX
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dcopy(int n,
+                 double const* vectorX,
+                 double *      vectorY)const{
+   int incrementX=1;
+   int incrementY=1;
+   this->Dcopy(n, vectorX, incrementX, vectorY, incrementY);
+}
+
+// vectorY = vectorX
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dcopy(int n,
+                 double const* vectorX, int incrementX,
+                 double*       vectorY, int incrementY) const{
+   double* x = const_cast<double*>(&vectorX[0]);
+   cblas_dcopy(n, x, incrementX, vectorY, incrementY);
+}
+
+// vectorY = alpha*vectorX + vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Daxpy(int n, double alpha,
+           double const* vectorX,
+           double*       vectorY) const{
+   int incrementX=1;
+   int incrementY=1;
+   this->Daxpy(n, alpha, vectorX, incrementX, vectorY, incrementY);
+}
+
+// vectorY = alpha*vectorX + vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Daxpy(int n, double alpha,
+           double const* vectorX, int incrementX,
+           double*       vectorY, int incrementY) const{
+   double* x = const_cast<double*>(&vectorX[0]);
+   cblas_daxpy(n, alpha, x, incrementX, vectorY, incrementY);
+}
+
+// returns vectorX^T*vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+double Blas::Ddot(int n,
+            double const* vectorX,
+            double const* vectorY) const{
+   int incrementX=1;
+   int incrementY=1;
+   return this->Ddot(n, vectorX, incrementX, vectorY, incrementY);
+}
+
+// returns vectorX^T*vectorY
+//    vectorX: n-vector
+//    vectorY: n-vector
+double Blas::Ddot(int n,
+            double const* vectorX, int incrementX,
+            double const* vectorY, int incrementY)const{
+   double* x=const_cast<double*>(vectorX),
+         * y=const_cast<double*>(vectorY);
+   return cblas_ddot(n, x, incrementX, y, incrementY);
+}
+
 // vectorY = matrixA*vectorX
 //    matrixA: m*n-matrix (matrixA[m][n] in row-major (C/C++ style))
 //    vectorX: n-vector
@@ -74,7 +138,7 @@ void Blas::Dgemv(int m, int n,
 }
 
 // vectorY = alpha*matrixA*vectorX + beta*vectorY
-//    matrixA: m*n-matrix (matrixA[m][n] in row-major (C/C++ style))
+//    matrixA: m*n-matrix 
 //    vectorX: n-vector
 //    vectorY: m-vector
 void Blas::Dgemv(bool isColumnMajorMatrixA,
@@ -101,6 +165,62 @@ void Blas::Dgemv(bool isColumnMajorMatrixA,
 
 }
 
+// vectorY = matrixA*vectorX
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part)
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dsymv(int n,
+           double const* const* matrixA,
+           double const* vectorX,
+           double*       vectorY) const{
+   bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major.
+   int incrementX=1, incrementY=1;
+   double alpha=1.0, beta=0.0;
+   this->Dsymv(n, alpha, matrixA, vectorX, incrementX, beta, vectorY, incrementY);
+}
+
+// vectorY = alpha*matrixA*vectorX + beta*vectorY
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part)
+//    vectorX: n-vector
+//    vectorY: n-vector
+void Blas::Dsymv(int n, double alpha,
+           double const* const* matrixA,
+           double const* vectorX, int incrementX,
+           double beta,
+           double*       vectorY, int incrementY) const{
+   double* a = const_cast<double*>(&matrixA[0][0]);
+   double* x = const_cast<double*>(&vectorX[0]);
+   CBLAS_UPLO uploA=CblasUpper;
+   int lda = n;
+   cblas_dsymv(CblasRowMajor, uploA, n, alpha, a, lda, x, incrementX, beta, vectorY, incrementY);
+}
+
+// matrixA = alpha*vectorX*vectorX^T + matrixA
+//    matrixA: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.)
+//    vectorX: n-matrix
+void Blas::Dsyr(int n, double alpha,
+          double const* vectorX,
+          double ** matrixA)const{
+   int incrementX=1;
+   this->Dsyr(n, alpha, vectorX, incrementX, matrixA);
+}
+
+void Blas::Dsyr(int n, double alpha,
+          double const* vectorX, int incrementX,
+          double ** matrixA)const{
+   double* a = &matrixA[0][0];
+   double* x = const_cast<double*>(&vectorX[0]);
+   CBLAS_UPLO uploA=CblasUpper;
+   int lda = n;
+   cblas_dsyr(CblasRowMajor, uploA, n, alpha, x, incrementX, a, lda);
+#pragma omp parallel for schedule(auto)
+   for(int i=0;i<n;i++){
+      for(int j=i+1;j<n;j++){
+         matrixA[j][i] = matrixA[i][j];
+      }
+   }
+}
+
 // matrixC = matrixA*matrixB
 //    matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style))
 //    matrixB: k*n-matrix (matrixB[k][n] in row-major (C/C++ style))
@@ -117,8 +237,8 @@ void Blas::Dgemm(int m, int n, int k,
 }
 
 // matrixC = alpha*matrixA*matrixB + beta*matrixC
-//    matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style))
-//    matrixB: k*n-matrix (matrixB[k][n] in row-major (C/C++ style))
+//    matrixA: m*k-matrix 
+//    matrixB: k*n-matrix 
 //    matrixC: m*n-matrix (matrixC[m][n] in row-major (C/C++ style))
 void Blas::Dgemm(bool isColumnMajorMatrixA, 
                  bool isColumnMajorMatrixB, 

--- a/src/wrappers/Lapack.cpp
+++ b/src/wrappers/Lapack.cpp
@@ -320,7 +320,16 @@ int Lapack::Dgetrs(double const* const* matrix, double** b, int size, int nrhs)
 // Argument "matrix" is sizeM*sizeN matrix.
 // Argument "matrix" will be LU-decomposed.
 int Lapack::Dgetrf(double** matrix, int sizeM, int sizeN) const{
-   int*    ipiv            = (int*)   mkl_malloc( sizeof(int)*2*sizeM,        16 );
+   int* ipiv = (int*)mkl_malloc( sizeof(int)*2*sizeM,        16 );
+   this->Dgetrf(matrix, ipiv, sizeM, sizeN);
+   mkl_free(ipiv);
+   int info = 0;
+   return info;
+}
+
+// Argument "matrix" is sizeM*sizeN matrix.
+// Argument "matrix" will be LU-decomposed.
+int Lapack::Dgetrf(double** matrix, int* ipiv, int sizeM, int sizeN) const{
    double* convertedMatrix = (double*)mkl_malloc( sizeof(double)*sizeM*sizeN, 16 );
    for(int i=0; i<sizeM; i++){
       for(int j=0; j<sizeN; j++){
@@ -334,7 +343,6 @@ int Lapack::Dgetrf(double** matrix, int sizeM, int sizeN) const{
       }
    }
    mkl_free(convertedMatrix);
-   mkl_free(ipiv);
    int info = 0;
    return info;
 }

--- a/src/wrappers/Lapack.h
+++ b/src/wrappers/Lapack.h
@@ -28,6 +28,7 @@ public:
    int Dsysv(double const* const* matrix, double* b, int size);
    int Dgetrs(double const* const* matrix, double** b, int size, int nrhs) const;
    int Dgetrf(double** matrix, int sizeM, int sizeN) const;
+   int Dgetrf(double** matrix, int* ipiv, int sizeM, int sizeN) const;
 private:
    Lapack();
    ~Lapack();

--- a/src/wrappers/Lapack_GNU.cpp
+++ b/src/wrappers/Lapack_GNU.cpp
@@ -319,7 +319,16 @@ int Lapack::Dgetrs(double const* const* matrix, double** b, int size, int nrhs)
 // Argument "matrix" is sizeM*sizeN matrix.
 // Argument "matrix" will be LU-decomposed.
 int Lapack::Dgetrf(double** matrix, int sizeM, int sizeN) const{
-   int*    ipiv            = (int*)   LAPACKE_malloc( sizeof(int)*2*sizeM );
+   int* ipiv = (int*)LAPACKE_malloc( sizeof(int)*2*sizeM );
+   this->Dgetrf(matrix, ipiv, sizeM, sizeN);
+   LAPACKE_free(ipiv);
+   int info = 0;
+   return info;
+}
+
+// Argument "matrix" is sizeM*sizeN matrix.
+// Argument "matrix" will be LU-decomposed.
+int Lapack::Dgetrf(double** matrix, int* ipiv, int sizeM, int sizeN) const{
    double* convertedMatrix = (double*)LAPACKE_malloc( sizeof(double)*sizeM*sizeN );
    for(int i=0; i<sizeM; i++){
       for(int j=0; j<sizeN; j++){
@@ -333,7 +342,6 @@ int Lapack::Dgetrf(double** matrix, int sizeM, int sizeN) const{
       }
    }
    LAPACKE_free(convertedMatrix);
-   LAPACKE_free(ipiv);
    int info = 0;
    return info;
 }

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -30,8 +30,10 @@
 #include"../base/PrintController.h"
 #include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
+#include"../wrappers/Blas.h"
 #include"../wrappers/Lapack.h"
 #include"../base/Enums.h"
+#include"../base/MathUtilities.h"
 #include"../base/MallocerFreer.h"
 #include"../base/EularAngle.h"
 #include"../base/Parameters.h"
@@ -696,6 +698,124 @@ void ZindoS::CalcDiatomicOverlapAOs2ndDerivativeInDiatomicFrame(double** diatomi
    diatomicOverlapAOs2ndDeri[px][px] *= this->overlapAOsCorrectionPi;
 }
 
+// calculate OverlapSingletSDs matrix between different electronic-structure, S^{SSD}_{ij}.
+// i and j are singlet SDs belonging to left and right hand side electronic-structures, respectively.
+// The index i=0 means the Hartree-Fock state.
+// This overlapSingletSDs are calculated from overlapMOs.
+// Note that rhs-electronic-structure is this electronic-structure  
+// and lhs-electronic-structure is another electronic-structure.
+void ZindoS::CalcOverlapSingletSDsWithAnotherElectronicStructure(double** overlapSingletSDs, 
+                                                                 double const* const* overlapMOs) const{
+   int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
+   double** tmpMatrix1=NULL;
+   double** tmpMatrix2=NULL;
+   double** tmpMatrix3=NULL;
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix1, numberOcc, numberOcc);
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix2, numberOcc, numberOcc);
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix3, numberOcc, numberOcc);
+   double sqrtGroundStateOverlap;
+   try{
+      // between ground state
+      for(int i=0; i<numberOcc; i++){
+         for(int j=0; j<numberOcc; j++){
+            tmpMatrix1[i][j] = overlapMOs[i][j];
+         }
+      }
+      sqrtGroundStateOverlap = GetDeterminant(tmpMatrix1, numberOcc);
+      overlapSingletSDs[0][0] = pow(sqrtGroundStateOverlap,2.0);
+
+      for(int k=0; k<this->matrixCISdimension; k++){
+         // single excitation from I-th (occupied)MO to A-th (virtual)MO
+         int moI = this->GetActiveOccIndex(*this->molecule, k);
+         int moA = this->GetActiveVirIndex(*this->molecule, k);
+         for(int l=0; l<this->matrixCISdimension; l++){
+            // single excitation from I-th (occupied)MO to A-th (virtual)MO
+            int moJ = this->GetActiveOccIndex(*this->molecule, l);
+            int moB = this->GetActiveVirIndex(*this->molecule, l);
+            for(int i=0; i<numberOcc; i++){
+               int destMO = i==moI ? moA : i;
+               for(int j=0; j<numberOcc; j++){
+                  int sourceMO = j==moJ ? moB : j;
+                  tmpMatrix1[i][j] = overlapMOs[destMO][j];
+                  tmpMatrix2[i][j] = overlapMOs[i][sourceMO];
+                  tmpMatrix3[i][j] = overlapMOs[destMO][sourceMO];
+               }
+            }
+            double det1 = GetDeterminant(tmpMatrix1, numberOcc);
+            double det2 = GetDeterminant(tmpMatrix2, numberOcc);
+            double det3 = GetDeterminant(tmpMatrix3, numberOcc);
+            // from ground state to singet SDs
+            overlapSingletSDs[k+1][0]   = sqrtGroundStateOverlap*det1;
+            // from singet SDs to ground state
+            overlapSingletSDs[0]  [l+1] = sqrtGroundStateOverlap*det2;
+            // from singet SDs to singlet SDs
+            overlapSingletSDs[k+1][l+1] = sqrtGroundStateOverlap*det3 + det1*det2;
+         }
+      }
+   }
+   catch(MolDSException ex){
+      MallocerFreer::GetInstance()->Free<double>(&tmpMatrix1, numberOcc, numberOcc);
+      MallocerFreer::GetInstance()->Free<double>(&tmpMatrix2, numberOcc, numberOcc);
+      MallocerFreer::GetInstance()->Free<double>(&tmpMatrix3, numberOcc, numberOcc);
+      throw ex;
+   }
+   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix1, numberOcc, numberOcc);
+   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix2, numberOcc, numberOcc);
+   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix3, numberOcc, numberOcc);
+}
+
+// calculate overlapESs (ES means eigenstate) matrix between different electronic-structure, S^{ES}_{ij}.
+// i and j are singlet SDs belonging to left and right hand side electronic-structures, respectively.
+// The index i=0 means the ground state.
+// This overlapESs is calculated from the overlapsingletSDs.
+// Note that rhs-electronic-structure is this electronic-structure  
+// and lhs-electronic-structure is another electronic-structure.
+void ZindoS::CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs, 
+                                                          double const* const* overlapSingletSDs,
+                                                          const MolDS_base::ElectronicStructure& lhsElectronicStructure) const{
+   const ElectronicStructure* rhsElectronicStructure = this;
+   double const* const* rhsMatrixCIS = this->matrixCIS;
+   double const* const* lhsMatrixCIS = lhsElectronicStructure.GetMatrixCIS();
+   int dimOverlapSingletSDs = this->matrixCISdimension + 1;
+   int dimOverlapESs = Parameters::GetInstance()->GetNumberElectronicStatesNASCO();
+   int groundstate = 0;
+   // extended CIS matrix includes groundstate althoug matrixCIS does not include groundstate.
+   double** lhsExtendedMatrixCIS=NULL;
+   double** rhsExtendedMatrixCIS=NULL;
+   double** tmpMatrix=NULL;
+   MallocerFreer::GetInstance()->Malloc<double>(&lhsExtendedMatrixCIS, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Malloc<double>(&rhsExtendedMatrixCIS, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,            dimOverlapSingletSDs, dimOverlapESs);
+   lhsExtendedMatrixCIS[groundstate][groundstate] = 1.0;
+   rhsExtendedMatrixCIS[groundstate][groundstate] = 1.0;
+   for(int i=1; i<dimOverlapESs; i++){
+      for(int j=1; j<dimOverlapSingletSDs; j++){
+         rhsExtendedMatrixCIS[i][j] = rhsMatrixCIS[i-1][j-1];
+         lhsExtendedMatrixCIS[i][j] = lhsMatrixCIS[i-1][j-1];
+      }
+   }
+   // calc. overlap between eigenstates
+   bool isColumnMajorOverlapSingletSDs = false;
+   bool isColumnMajorRhsMatrixCIS = true;
+   double alpha=1.0;
+   double beta=0.0;
+   MolDS_wrappers::Blas::GetInstance()->Dgemm(isColumnMajorOverlapSingletSDs,
+                                              isColumnMajorRhsMatrixCIS,
+                                              dimOverlapSingletSDs, dimOverlapESs, dimOverlapSingletSDs,
+                                              alpha,
+                                              overlapSingletSDs,
+                                              rhsExtendedMatrixCIS,
+                                              beta,
+                                              tmpMatrix);
+   MolDS_wrappers::Blas::GetInstance()->Dgemm(dimOverlapESs, dimOverlapESs, dimOverlapSingletSDs,
+                                              lhsExtendedMatrixCIS,
+                                              tmpMatrix,
+                                              overlapESs);
+   MallocerFreer::GetInstance()->Free<double>(&lhsExtendedMatrixCIS, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Free<double>(&rhsExtendedMatrixCIS, dimOverlapSingletSDs, dimOverlapSingletSDs);
+   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,            dimOverlapSingletSDs, dimOverlapESs);
+}
+
 // The order of mol, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
 double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL, 
                                            const Molecule& molecule, 

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -29,6 +29,11 @@ public:
    virtual ~ZindoS();
    void DoCIS();
    void OutputCISResults() const;
+   void CalcOverlapSingletSDsWithAnotherElectronicStructure(double** overlapSingletSDs, 
+                                                            double const* const* overlapMOs) const;
+   void CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs, 
+                                                     double const* const* overlapSingletSDs,
+                                                     const MolDS_base::ElectronicStructure& lhsElectronicStructure) const;
 protected:
    std::string errorMessageDavidsonNotConverged;
    std::string errorMessageCalcCISMatrix;