Commit: ce896c8ce3c7ae8cbfadc32e4101add78c19e25f - MolDS (git)

Commit MetaInfo

• Révision: ce896c8ce3c7ae8cbfadc32e4101add78c19e25f (tree)
• l'heure: 2011-07-08 19:20:41
• Auteur: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Message de Log

CIS-Davidson is implemented. This is implemented with also open MP.

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

Change Summary

• modified: src/base/InputParser.h
• modified: src/base/Parameters.h
• modified: src/input.in
• modified: src/mkl_wrapper/LapackWrapper.h
• modified: src/zindo/ZindoS.h

Modification

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -48,6 +48,9 @@ private:
    string messageCisNumberActiveVir;
    string messageCisNumberExcitedStates;
    string messageCisDavidson;
+   string messageCisNormTolerance;
+   string messageCisMaxIterations;
+   string messageCisMaxDimensions;
    string stringYES;
    string stringNO;
    string stringSpace;
@@ -93,6 +96,9 @@ private:
    string stringCISActiveVir;
    string stringCISNStates;
    string stringCISDavidson;
+   string stringCISMaxIter;
+   string stringCISMaxDimensions;
+   string stringCISNormTolerance;
    void CalcMolecularBasics(Molecule* molecule);
    void CalcCisCondition(Molecule* molecule);
    void OutputMolecularBasics(Molecule* molecule);
@@ -124,6 +130,9 @@ InputParser::InputParser(){
    this->messageCisNumberActiveVir = "\t\tNumber of active Vir.: ";
    this->messageCisNumberExcitedStates = "\t\tNumber of excited states: ";
    this->messageCisDavidson = "\t\tCIS-Davidson: ";
+   this->messageCisNormTolerance = "\t\tNorm tolerance for the residual of the Davidson: ";
+   this->messageCisMaxIterations = "\t\tMax iterations for the Davidson: ";
+   this->messageCisMaxDimensions = "\t\tMax dimensions for the Davidson: ";
    this->stringYES = "yes";
    this->stringNO = "no";
    this->stringSpace = " ";
@@ -169,6 +178,9 @@ InputParser::InputParser(){
    this->stringCISActiveVir = "active_vir";
    this->stringCISNStates = "nstates";
    this->stringCISDavidson = "davidson";
+   this->stringCISMaxIter = "max_iter";
+   this->stringCISMaxDimensions = "max_dim";
+   this->stringCISNormTolerance = "norm_tol";
 }
 
 InputParser::~InputParser(){
@@ -426,12 +438,29 @@ void InputParser::Parse(Molecule* molecule){
                }
                j++;
             }
+            // max iterations for the Davidson roop
+            if(inputTerms[j].compare(this->stringCISMaxIter) == 0){
+               int maxIter = atoi(inputTerms[j+1].c_str());
+               Parameters::GetInstance()->SetMaxIterationsCIS(maxIter);
+               j++;
+            }
+            // max dimensions for the Davidson expansion
+            if(inputTerms[j].compare(this->stringCISMaxDimensions) == 0){
+               int maxDim = atoi(inputTerms[j+1].c_str());
+               Parameters::GetInstance()->SetMaxDimensionsCIS(maxDim);
+               j++;
+            }
+            // nolm tolerance for the norm of the resiudal vectors of the Davidson.
+            if(inputTerms[j].compare(this->stringCISNormTolerance) == 0){
+               double normTol = atof(inputTerms[j+1].c_str());
+               Parameters::GetInstance()->SetNormToleranceCIS(normTol);
+               j++;
+            }
             j++;   
          }
          i = j;
       }
 
-
       // theory
       if(inputTerms[i].compare(this->stringTheory) == 0){
          int j=i+1;
@@ -524,14 +553,17 @@ void InputParser::CalcCisCondition(Molecule* molecule){
    }   
 
    // check the number of calculated excited states.
-   int numberExcitedStates = Parameters::GetInstance()->GetActiveOccCIS() 
-                            *Parameters::GetInstance()->GetActiveVirCIS();
+   int numberSlaterDeterminants = Parameters::GetInstance()->GetActiveOccCIS() 
+                                 *Parameters::GetInstance()->GetActiveVirCIS();
    if(!Parameters::GetInstance()->GetIsDavidsonCIS()){
-      Parameters::GetInstance()->SetNumberExcitedStatesCIS(numberExcitedStates);
+      Parameters::GetInstance()->SetNumberExcitedStatesCIS(numberSlaterDeterminants);
    }
    else{
-      if(numberExcitedStates < Parameters::GetInstance()->GetNumberExcitedStatesCIS()){
-         Parameters::GetInstance()->SetNumberExcitedStatesCIS(numberExcitedStates);
+      if(numberSlaterDeterminants < Parameters::GetInstance()->GetNumberExcitedStatesCIS()){
+         Parameters::GetInstance()->SetNumberExcitedStatesCIS(numberSlaterDeterminants);
+      }
+      if(numberSlaterDeterminants < Parameters::GetInstance()->GetMaxDimensionsCIS()){
+         Parameters::GetInstance()->SetMaxDimensionsCIS(numberSlaterDeterminants);
       }
    }
    
@@ -568,6 +600,9 @@ void InputParser::OutputCisConditions(){
    printf("%s",this->messageCisDavidson.c_str());
    if(Parameters::GetInstance()->GetIsDavidsonCIS()){
       printf("%s\n",this->stringYES.c_str());
+      printf("%s%d\n",this->messageCisMaxIterations.c_str(),Parameters::GetInstance()->GetMaxIterationsCIS());
+      printf("%s%d\n",this->messageCisMaxDimensions.c_str(),Parameters::GetInstance()->GetMaxDimensionsCIS());
+      printf("%s%e\n",this->messageCisNormTolerance.c_str(),Parameters::GetInstance()->GetNormToleranceCIS());
    }
    else{
       printf("%s\n",this->stringNO.c_str());

--- a/src/base/Parameters.h
+++ b/src/base/Parameters.h
@@ -59,6 +59,12 @@ public:
    void SetNumberExcitedStatesCIS(int nStates);
    bool GetIsDavidsonCIS();
    void SetIsDavidsonCIS(bool isDavidsonCIS);
+   int GetMaxIterationsCIS();
+   void SetMaxIterationsCIS(int maxIterationsCIS);
+   int GetMaxDimensionsCIS();
+   void SetMaxDimensionsCIS(int maxDimensionsCIS);
+   double GetNormToleranceCIS();
+   void SetNormToleranceCIS(double normToleranceCIS);
 private:
    static Parameters* parameters;
    Parameters();
@@ -91,6 +97,9 @@ private:
    int activeOccCIS;
    int activeVirCIS;
    int numberExcitedStatesCIS;
+   int maxIterationsCIS;
+   int maxDimensionsCIS;
+   double normToleranceCIS;
    bool isDavidsonCIS;
 };
 Parameters* Parameters::parameters = NULL;
@@ -155,6 +164,9 @@ void Parameters::SetDefaultValues(){
    this->activeVirCIS = 10;
    this->numberExcitedStatesCIS = 5;
    this->isDavidsonCIS = true;
+   this->maxIterationsCIS = 100;
+   this->maxDimensionsCIS = 100;
+   this->normToleranceCIS = pow(10.0, -6.0);
 }
 
 double Parameters::GetThresholdSCF(){
@@ -354,6 +366,30 @@ void Parameters::SetIsDavidsonCIS(bool isDavidsonCIS){
    this->isDavidsonCIS = isDavidsonCIS;
 }
 
+int Parameters::GetMaxIterationsCIS(){
+   return this->maxIterationsCIS;
+}
+
+void Parameters::SetMaxIterationsCIS(int maxIterationsCIS){
+   this->maxIterationsCIS = maxIterationsCIS;
+}
+
+int Parameters::GetMaxDimensionsCIS(){
+   return this->maxDimensionsCIS;
+}
+
+void Parameters::SetMaxDimensionsCIS(int maxDimensionsCIS){
+   this->maxDimensionsCIS = maxDimensionsCIS;
+}
+
+double Parameters::GetNormToleranceCIS(){
+   return this->normToleranceCIS;
+}
+
+void Parameters::SetNormToleranceCIS(double normToleranceCIS){
+   this->normToleranceCIS = normToleranceCIS;
+}
+
 }
 #endif
 

--- a/src/input.in
+++ b/src/input.in
@@ -22,7 +22,10 @@ CIS
    davidson yes
    active_occ 100
    active_vir 100
-   nstates 10
+   nstates 11
+   max_iter 200
+   max_dim 100
+   norm_tol 0.000001
 CIS_END
 
 
@@ -80,13 +83,13 @@ CIS_END
 //GEOMETRY_END
 
 //metane
-GEOMETRY
-   C  0.647389    0.820131    0.000000
-   H  1.004043   -0.188679    0.000000
-   H  1.004062    1.324529    0.873652
-   H  1.004062    1.324529   -0.873652
-   H  -0.422611    0.820144    0.000000
-GEOMETRY_END
+//GEOMETRY
+//   C  0.647389    0.820131    0.000000
+//   H  1.004043   -0.188679    0.000000
+//   H  1.004062    1.324529    0.873652
+//   H  1.004062    1.324529   -0.873652
+//   H  -0.422611    0.820144    0.000000
+//GEOMETRY_END
 
 
 // c2
@@ -131,12 +134,12 @@ GEOMETRY_END
 //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
+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

--- a/src/mkl_wrapper/LapackWrapper.h
+++ b/src/mkl_wrapper/LapackWrapper.h
@@ -37,7 +37,7 @@ LapackWrapper* LapackWrapper::lapackWrapper = NULL;
 LapackWrapper::LapackWrapper(){
    this->calculatedDsysvBlockSize = false;
    this->dsysvBlockSize = 64;
-   this->errorMessageDsyevdInfo = "Error in mkl_wrapper::LapackWrapper::Dsyevd: info != 0\n";
+   this->errorMessageDsyevdInfo = "Error in mkl_wrapper::LapackWrapper::Dsyevd: info != 0: info = ";
    this->errorMessageDsyevdSize = "Error in mkl_wrapper::LapackWrapper::Dsyevd: size of matirx < 1\n";
    this->errorMessageDsysvInfo = "Error in mkl_wrapper::LapackWrapper::Dsysv: info != 0\n";
    this->errorMessageDsysvSize = "Error in mkl_wrapper::LapackWrapper::Dsysv: size of matirx < 1\n";
@@ -160,6 +160,7 @@ int LapackWrapper::Dsyevd(double** matrix, double* eigenValues, int size, bool c
    if(info != 0){
       stringstream ss;
       ss << errorMessageDsyevdInfo;
+      ss << info << endl;
       throw MolDSException(ss.str());
    }
    return info;

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -48,12 +48,46 @@ private:
                          Atom* atom); // Apendix in [BZ_1979]
    double GetNishimotoMatagaTwoEleInt(Atom* atomA, OrbitalType orbitalA, 
                                            Atom* atomB, OrbitalType orbitalB); // ref. [MN_1957] and (5a) in [AEZ_1986]
+   struct MoEnergy{
+      double energy;
+      int occIndex;
+      int virIndex;
+      int slaterIndex;
+   };
+   struct LessMoEnergy { bool operator()(const MoEnergy& rLeft, const MoEnergy& rRight) 
+      const { return rLeft.energy < rRight.energy; } };
    void DoesCISDirect();
    void DoesCISDavidson();
+   void CalcRitzVector(double* ritzVector, double** expansionVectors, double** interactionMatrix, 
+                       int interactionMatrixDimension, int ritzVectorIndex);
+   void CalcResidualVectorAndNorm(double* residualVector, double* norm, double* ritzVector, 
+                                    double* interactionEigenEnergies, int residualVectorIndex);
+   void SortSingleExcitationSlaterDeterminants(vector<MoEnergy>* moEnergies);
+   void UpdateExpansionVectors(double** expansionVectors, double* interactionEigenEnergies, double* residualVector,
+                               int interactionMatrixDimension, int* notConvergedStates, int residualVectorIndex);
+   void CalcInteractionMatrix(double** interactionMatrix, double** expansionVectors, int interactionMatrixDimension);
+   void FreeDavidsonCISTemporaryMtrices(double*** expansionVectors, double** residualVector, double** ritzVector);
+   void FreeDavidsonRoopCISTemporaryMtrices(double*** interactionMatrix, 
+                                        double interactionMatrixDimension, 
+                                        double** interactionEigenEnergies);
    void CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir);
    string errorMessageNishimotoMataga;
+   string errorMessageDavidsonNotConverged;
+   string errorMessageDavidsonMaxIter;
+   string errorMessageDavidsonMaxDim;
    string messageStartCIS;
    string messageDoneCIS;
+   string messageStartDirectCIS;
+   string messageDoneDirectCIS;
+   string messageStartDavidsonCIS;
+   string messageDoneDavidsonCIS;
+   string messageNumIterCIS;
+   string messageResidualNorm;
+   string messageDavidsonConverge;
+   string messageDavidsonReachCISMatrix;
+   string messageDavidsonGoToDirect;
+   string messageExcitedStatesEnergies;
+   string messageExcitedStatesEnergiesTitle;
 };
 
 ZindoS::ZindoS() : MolDS_cndo::Cndo2(){
@@ -92,11 +126,25 @@ void ZindoS::SetMessages(){
    this->errorMessageNishimotoMataga = "Error in base_zindo::ZindoS::GetNishimotoMatagaTwoEleInt: Invalid orbitalType.\n";
    this->errorMessageMolecularIntegralElement
       = "Error in zindo::ZindoS::GetMolecularIntegralElement: Non available orbital is contained.\n";
+   this->errorMessageDavidsonNotConverged =  "Error in zindo::ZindoS::DoesCISDavidson: Davidson did not met convergence criterion. \n";
+   this->errorMessageDavidsonMaxIter = "Davidson roop reaches max_iter=";
+   this->errorMessageDavidsonMaxDim = "Dimension of the expansion vectors reaches max_dim=";
    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";
    this->messageStartCIS = "**********  START: ZINDO/S-CIS  **********\n";
    this->messageDoneCIS = "**********  DONE: ZINDO/S-CIS  **********\n\n\n";
+   this->messageStartDirectCIS = "\t======  START: Direct-CIS  =====\n\n";
+   this->messageDoneDirectCIS =  "\t======  DONE: Direct-CIS  =====\n\n\n";
+   this->messageStartDavidsonCIS = "\t======  START: Davidson-CIS  =====\n";
+   this->messageDoneDavidsonCIS =  "\t======  DONE: Davidson-CIS  =====\n\n\n";
+   this->messageNumIterCIS = "\tDavidson iter=";
+   this->messageResidualNorm = "-th excited: norm of the residual = ";
+   this->messageDavidsonConverge = "\n\n\t\tDavidson for ZINDO/S-CIS met convergence criterion(^^b\n\n\n";
+   this->messageDavidsonReachCISMatrix = "\n\t\tDimension of the expansion vectors reaches to the dimension of the CIS-matrix.\n";
+   this->messageDavidsonGoToDirect = "\t\tHence, we go to the Direct-CIS.\n\n";
+   this->messageExcitedStatesEnergies = "\tExcitation energies:\n";
+   this->messageExcitedStatesEnergiesTitle = "\t\t| i-th | e[a.u.] | e[eV] | \n";
 }
 
 void ZindoS::SetEnableAtomTypes(){
@@ -621,6 +669,33 @@ double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
 void ZindoS::DoesCIS(){
    cout << this->messageStartCIS;
 
+   int numberOcc = Parameters::GetInstance()->GetActiveOccCIS();
+   int numberVir = Parameters::GetInstance()->GetActiveVirCIS();
+
+   // malloc or initialize  CIS matrix
+   if(this->matrixCIS == NULL){
+      this->matrixCISdimension = numberOcc*numberVir;
+      this->matrixCIS = MallocerFreer::GetInstance()->MallocDoubleMatrix2d(this->matrixCISdimension, 
+                                                                           this->matrixCISdimension);
+   }
+   else{
+      MallocerFreer::GetInstance()->InitializeDoubleMatrix2d(this->matrixCIS, 
+                                                             this->matrixCISdimension, 
+                                                             this->matrixCISdimension);
+   }
+
+   // malloc or initialize CIS eigen vector
+   if(this->excitedEnergies == NULL){
+      this->excitedEnergies = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(this->matrixCISdimension);
+   }
+   else{
+      MallocerFreer::GetInstance()->InitializeDoubleMatrix1d(this->excitedEnergies, this->matrixCISdimension);
+   }
+
+   // calculate CIS matrix
+   this->CalcCISMatrix(matrixCIS, numberOcc, numberVir);
+
+   // calculate excited energies
    if(Parameters::GetInstance()->GetIsDavidsonCIS()){
       this->DoesCISDavidson();
    }
@@ -629,55 +704,325 @@ void ZindoS::DoesCIS(){
    }
 
    // output eigen energies
+   cout << this->messageExcitedStatesEnergies;
+   cout << this->messageExcitedStatesEnergiesTitle;
+   double eV2AU = Parameters::GetInstance()->GetEV2AU();
    for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
-      printf("%d-th excited energy: %e\n",k+1, this->excitedEnergies[k]);
+      printf("\t\t %d\t%e\t%e\n",k+1, this->excitedEnergies[k], this->excitedEnergies[k]/eV2AU);
    }
    cout << endl;
    cout << this->messageDoneCIS;
 }
 
-void ZindoS::DoesCISDavidson(){
-   // ToDo: CIS-Davidson
-   stringstream ss;
-   ss << "CIS-Davidson (ZindoS::DoesCISDavidson()) is not implemented yet.\n\n";
-   throw MolDSException(ss.str());
+void ZindoS::SortSingleExcitationSlaterDeterminants(vector<MoEnergy>* moEnergies){
+
+   int numberOcc = Parameters::GetInstance()->GetActiveOccCIS();
+   int numberVir = Parameters::GetInstance()->GetActiveVirCIS();
+   for(int k=0; k<numberOcc*numberVir; k++){
+      // single excitation from I-th (occupied)MO to A-th (virtual)MO
+      int moI = this->molecule->GetTotalNumberValenceElectrons()/2 - (k/numberVir) -1;
+      int moA = this->molecule->GetTotalNumberValenceElectrons()/2 + (k%numberVir);
+
+      MoEnergy moEnergy = {this->energiesMO[moA] - this->energiesMO[moI], moI, moA, k};
+      moEnergies->push_back(moEnergy);
+   }
+   sort(moEnergies->begin(), moEnergies->end(), LessMoEnergy());
 }
 
-void ZindoS::DoesCISDirect(){
+// This method is used for Davidson
+void ZindoS::CalcRitzVector(double* ritzVector, double** expansionVectors, double** interactionMatrix, 
+                            int interactionMatrixDimension, int ritzVectorIndex){
+
+   for(int j=0; j<this->matrixCISdimension; j++){
+      ritzVector[j] = 0.0;
+      for(int k=0; k<interactionMatrixDimension; k++){
+         ritzVector[j] += expansionVectors[j][k]*interactionMatrix[ritzVectorIndex][k];
+      }
+   }
+}
+
+// This method is used for Davidson
+void ZindoS::CalcResidualVectorAndNorm(double* residualVector, double* norm, double* ritzVector, 
+                                       double* interactionEigenEnergies, int residualVectorIndex){
+
+   double sqNorm = 0.0;
+   for(int j=0; j<this->matrixCISdimension; j++){
+      residualVector[j] = interactionEigenEnergies[residualVectorIndex] * ritzVector[j];
+      for(int k=0; k<this->matrixCISdimension; k++){
+         double value = j<=k ? this->matrixCIS[j][k] : this->matrixCIS[k][j];
+         residualVector[j] -= value*ritzVector[k];
+      }
+      sqNorm += pow(residualVector[j],2.0);
+   }
+   *norm = sqrt(sqNorm);
+}
+
+// This method is used for Davidson
+void ZindoS::UpdateExpansionVectors(double** expansionVectors, double* interactionEigenEnergies, double* residualVector,
+                                    int interactionMatrixDimension, int* notConvergedStates, int residualVectorIndex){
+
+   double newExpansionVector[this->matrixCISdimension];
+
+   // calculate new expansion vector from residual vector
+   for(int j=0; j<this->matrixCISdimension; j++){
+      double temp = interactionEigenEnergies[residualVectorIndex]-this->matrixCIS[j][j];
+      if(temp == 0.0){
+         // prevent dividing by 0.
+         temp = pow(10,-100);
+      }
+      newExpansionVector[j]=pow(temp, -1.0)*residualVector[j];
+   }
+
+   // orthonormalize old expansion vectors and new expansion vector
+   for(int k=0; k<interactionMatrixDimension+*notConvergedStates; k++){
+      double overlap=0.0;
+      for(int j=0; j<this->matrixCISdimension; j++){
+         overlap += expansionVectors[j][k] * newExpansionVector[j];
+      }
+      for(int j=0; j<this->matrixCISdimension; j++){
+         newExpansionVector[j] -= overlap*expansionVectors[j][k];
+      }
+   }
+
+   // add new expansion vector to old expansion vectors
+   double sqNormNewExpVect = 0.0;
+   for(int j=0; j<this->matrixCISdimension; j++){
+      sqNormNewExpVect += pow(newExpansionVector[j],2.0);
+   }
+   double normNewExpVect = sqrt(sqNormNewExpVect);
+   for(int j=0; j<this->matrixCISdimension; j++){
+      expansionVectors[j][interactionMatrixDimension+*notConvergedStates] 
+            = pow(normNewExpVect,-1.0)*newExpansionVector[j];
+   }
+   *notConvergedStates += 1;
+}
+
+// This method is used for Davidson
+void ZindoS::CalcInteractionMatrix(double** interactionMatrix, double** expansionVectors, int interactionMatrixDimension){
+
+   for(int i=0; i<interactionMatrixDimension; i++){
+      for(int j=i; j<interactionMatrixDimension; j++){
+         for(int a=0; a<this->matrixCISdimension; a++){
+            for(int b=0; b<this->matrixCISdimension; b++){
+               double value = a<=b ? this->matrixCIS[a][b] : this->matrixCIS[b][a];
+               interactionMatrix[i][j] += expansionVectors[a][i] 
+                                          *value
+                                          *expansionVectors[b][j];
+            }
+         }
+      }
+   }
+}
+
+void ZindoS::DoesCISDavidson(){
+   cout << this->messageStartDavidsonCIS;
 
    int numberOcc = Parameters::GetInstance()->GetActiveOccCIS();
    int numberVir = Parameters::GetInstance()->GetActiveVirCIS();
    int numberExcitedStates = Parameters::GetInstance()->GetNumberExcitedStatesCIS();
+   int maxIter = Parameters::GetInstance()->GetMaxIterationsCIS();
+   int maxDim  = Parameters::GetInstance()->GetMaxDimensionsCIS();
+   double normTol = Parameters::GetInstance()->GetNormToleranceCIS();
+   bool convergeExcitedStates[numberExcitedStates];
+   int interactionMatrixDimension;
+   bool reachMaxDim;
+   bool allConverged;
+   int notConvergedStates;
+   bool goToDirectCIS;
+   double** expansionVectors = NULL;
+   double*  ritzVector = NULL;
+   double*  residualVector = NULL;
+
+   // malloc
+   expansionVectors = MallocerFreer::GetInstance()->MallocDoubleMatrix2d(this->matrixCISdimension, maxDim);
+   ritzVector = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(this->matrixCISdimension);
+   residualVector = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(this->matrixCISdimension);
+
+   try{
+      // sort single excitation slater determinants
+      vector<MoEnergy> moEnergies;
+      this->SortSingleExcitationSlaterDeterminants(&moEnergies);
+
+      // set initial expansion vectors and initial conveged vectors
+      for(int k=0; k<numberExcitedStates; k++){
+         expansionVectors[moEnergies[k].slaterIndex][k] = 1.0;
+         convergeExcitedStates[k] = false;
+      }
 
-   // malloc CIS matrix
-   if(this->matrixCIS == NULL){
-      this->matrixCISdimension = numberExcitedStates;
-      this->matrixCIS = MallocerFreer::GetInstance()->MallocDoubleMatrix2d(this->matrixCISdimension, 
-                                                                           this->matrixCISdimension);
+      interactionMatrixDimension = 0;
+      reachMaxDim = false;
+      goToDirectCIS = false;
+      // Davidson roop
+      for(int k=0; k<maxIter; k++){
+         cout << messageNumIterCIS << k << endl;
+
+         // calculate dimension of the interaction matrix (= number of the expansion vectors).
+         for(int i=0; i<numberExcitedStates; i++){
+            if(!convergeExcitedStates[i]){
+               interactionMatrixDimension += 1;
+            }
+         }
+
+         // malloc interaction matrix and etc.
+         double** interactionMatrix = NULL;
+         double* interactionEigenEnergies = NULL;
+         interactionMatrix = MallocerFreer::GetInstance()->MallocDoubleMatrix2d
+                              (interactionMatrixDimension, interactionMatrixDimension);
+         interactionEigenEnergies = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(interactionMatrixDimension);
+
+         try{
+
+            // calculate interaction matrix
+            this->CalcInteractionMatrix(interactionMatrix, expansionVectors, interactionMatrixDimension);
+
+            // diagonalize interaction matrix
+            bool calcEigenVectors = true;
+            MolDS_mkl_wrapper::LapackWrapper::GetInstance()->Dsyevd(interactionMatrix,
+                                                                    interactionEigenEnergies, 
+                                                                    interactionMatrixDimension, 
+                                                                    calcEigenVectors);
+
+            // check convergence of all excited states
+            notConvergedStates=0;
+            allConverged = true;
+            for(int i=0; i<numberExcitedStates; i++){
+       
+               // calculate i-th ritz vector
+               this->CalcRitzVector(ritzVector, expansionVectors, interactionMatrix, interactionMatrixDimension, i);
+
+               // calculate i-th residual vector and the norm of the residual vector
+               double norm = 0.0;
+               this->CalcResidualVectorAndNorm(residualVector, &norm, ritzVector, interactionEigenEnergies, i);
+
+               // output norm of residual vector
+               cout << "\t  " << i+1 << this->messageResidualNorm << norm << endl;
+               if(i == numberExcitedStates-1){
+                  cout << endl;
+               }
+        
+               // check tolerance for the norm of the residual vector.
+               if(norm < normTol){
+                  convergeExcitedStates[i] = true;
+               }
+               else{
+                  convergeExcitedStates[i] = false;
+                  allConverged = false;
+                  if(interactionMatrixDimension+notConvergedStates == maxDim && maxDim !=this->matrixCISdimension){
+                     reachMaxDim = true;
+                     break;
+                  }
+                  else if(interactionMatrixDimension+notConvergedStates == this->matrixCISdimension){
+                     goToDirectCIS = true;
+                     break;
+                  }
+            
+                  // update expansion vectors
+                  this->UpdateExpansionVectors(expansionVectors, interactionEigenEnergies, residualVector,
+                                                interactionMatrixDimension, &notConvergedStates, i);
+
+               }
+            } 
+
+            if(allConverged){
+               // copy to cis eigen vector and value
+               for(int i=0; i<numberExcitedStates; i++){
+                  this->excitedEnergies[i] = interactionEigenEnergies[i];
+                  this->CalcRitzVector(ritzVector, expansionVectors, interactionMatrix, interactionMatrixDimension, i);
+                  for(int j=0; j<this->matrixCISdimension; j++){
+                     this->matrixCIS[i][j] = ritzVector[j];
+                  }
+               }
+            }
+         }
+         catch(MolDSException ex){
+            this->FreeDavidsonRoopCISTemporaryMtrices(&interactionMatrix, 
+                                                      interactionMatrixDimension, 
+                                                      &interactionEigenEnergies);
+            throw ex;
+         }
+         this->FreeDavidsonRoopCISTemporaryMtrices(&interactionMatrix, 
+                                                   interactionMatrixDimension, 
+                                                   &interactionEigenEnergies);
+
+         // stop the Davidson roop
+         if(allConverged){
+            cout << this->messageDavidsonConverge;
+            break;
+         }
+         else if(!allConverged && goToDirectCIS){
+            cout << this->messageDavidsonReachCISMatrix;
+            cout << this->messageDavidsonGoToDirect;
+            break;
+         }
+         else if(!allConverged && reachMaxDim){
+            stringstream ss;
+            ss << endl;
+            ss << this->errorMessageDavidsonNotConverged;
+            ss << this->errorMessageDavidsonMaxDim << maxDim << endl;
+            throw MolDSException(ss.str());
+         }
+         else if(!allConverged && k==maxIter-1){
+            stringstream ss;
+            ss << this->errorMessageDavidsonNotConverged;
+            ss << this->errorMessageDavidsonMaxIter << maxIter << endl;
+            throw MolDSException(ss.str());
+         }
+
+      }// end Davidson roop
    }
-   else{
-      MallocerFreer::GetInstance()->InitializeDoubleMatrix2d(this->matrixCIS, 
-                                                             this->matrixCISdimension, 
-                                                             this->matrixCISdimension);
+   catch(MolDSException ex){
+      this->FreeDavidsonCISTemporaryMtrices(&expansionVectors, &residualVector, &ritzVector);
+      throw ex;
    }
+   this->FreeDavidsonCISTemporaryMtrices(&expansionVectors, &residualVector, &ritzVector);
 
-   // malloc CIS eigen vector
-   if(this->excitedEnergies == NULL){
-      this->excitedEnergies = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(this->matrixCISdimension);
+   cout << this->messageDoneDavidsonCIS;
+   // change algorithm from Davidso to direct
+   if(goToDirectCIS){
+      this->DoesCISDirect();
    }
-   else{
-      MallocerFreer::GetInstance()->InitializeDoubleMatrix1d(this->excitedEnergies, this->matrixCISdimension);
+
+}
+
+void ZindoS::FreeDavidsonCISTemporaryMtrices(double*** expansionVectors, double** residualVector, double** ritzVector){
+
+   if(*expansionVectors != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix2d(expansionVectors, this->matrixCISdimension);
+      //cout << "expansion vectors deleted\n";
+   }
+   if(*residualVector != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(residualVector);
+      //cout << "residual vector deleted\n";
+   }
+   if(*ritzVector != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(ritzVector);
+      //cout << "ritz vector deleted\n";
    }
 
-   // calc CIS matrix
-   this->CalcCISMatrix(matrixCIS, numberOcc, numberVir);
+}
 
+void ZindoS::FreeDavidsonRoopCISTemporaryMtrices(double*** interactionMatrix, 
+                                             double interactionMatrixDimension, 
+                                             double** interactionEigenEnergies){
+   if(*interactionMatrix != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix2d(interactionMatrix, interactionMatrixDimension);
+      //cout << "interactionMatrix deleted\n";
+   }
+   if(*interactionEigenEnergies != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(interactionEigenEnergies);
+      //cout << "interactionEigenEnergies deleted\n";
+   }
+
+}
+void ZindoS::DoesCISDirect(){
+
+   cout << this->messageStartDirectCIS;
    bool calcEigenVectors = true;
    MolDS_mkl_wrapper::LapackWrapper::GetInstance()->Dsyevd(this->matrixCIS,
                                                            this->excitedEnergies, 
                                                            this->matrixCISdimension, 
                                                            calcEigenVectors);
-
+   cout << this->messageDoneDirectCIS;
 }
 
 void ZindoS::CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir){