提交: 0c1e048209b1eaecdb0fd329442b2c8a8bff4461 - MolDS (git)

Commit MetaInfo

• 修订版: 0c1e048209b1eaecdb0fd329442b2c8a8bff4461 (tree)
• 时间: 2011-02-02 11:57:29
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

damp and diis are added.

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

更改概述

• modified: doc/README.txt
• modified: doc/refferences.txt
• modified: src/base/InputParser.h
• modified: src/base/Parameters.h
• modified: src/cndo/Cndo2.h
• modified: src/input.in
• modified: src/mkl_wrapper/LapackWrapper.h

差异

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -20,23 +20,27 @@ SCF:
          THEORY_END
    
    -options
-    "max_iter", "rms_density", "damping_thresh", and "damping_weight" are prepared as options.
-    option is "origin" only for setting the origin of the inertia tensor.
-    options are written in inertia-directive in angstrom unit.
-    Center of mass is used as origin when the "origin" is not set.
+    "max_iter", "rms_density", "damping_thresh", "damping_weight", 
+    "diis_num_error_vect", "diis_start_error", and "diis_end_error" are prepared as options.
 
     Default value of "max_iter" is 100.
     Default value of "rms_density" is 10**(-8.0).
-    Default value of "dampingThresh" is 1.
-    Default value of "dampingWeight" is 0.8.
+    Default value of "damping_thresh" is 1.
+    Default value of "damping_weight" is 0.8.
+    Default value of "diis_num_error_vect" is 5.
+    Default value of "diis_start_error" is 0.01.
+    Default value of "diis_end_error" is 10**(-8.0).
 
     E.g.
-         INERTIA
+         SCF
             max_iter 200
             rms_density 0.00000001
             damping_thresh 0.1
             damping_weight 0.7
-         INERTIA_END
+            diis_num_error_vect 6
+            diis_start_error 0.01
+            diis_end_error 0.00000001
+         SCF_END
 
 Principal Axes (Diagonalizing the inertia tensor):
    Write "principal_axes" in theory-directive.

--- a/doc/refferences.txt
+++ b/doc/refferences.txt
@@ -7,6 +7,7 @@
 [RZ_1973] J. Ridley and M. C. Zerner, Theort. Chim. Acta (Berl.) 32, 111 (1973)
 [RZ_1976] J. E. Ridley and M. C. Zerner, Theort. Chim. Acta (Berl.) 42, 223 (1976)
 [BZ_1979] A. D. Bacon and M. C. Zerner, Theort. Chim. Acta (Berl.) 53, 21 (1979)
+[P_1980] P. Pulay, Chem. Phys. Lett. 73, 393 (1980)
 [HKLWNZ_1982] Z. S. Herman, R. F. Kirachner, G. H. Loew, U. T. M.-WesterHoff, A. Nazzal, and M. C. Zerner, Inorg. Chem., 21, 46 (1982)
 [AEZ_1986] W. P. Anderson, W. D. Edwards, and M. C. Zerner, Inorg. Chem. 25, 2728 (1986)
 [BFB_1997] M. A. Blanco, M. Fl{\'o}rez, and M. Bermejo, J. Mol. Struct. 419, 19 (1997)

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -40,6 +40,9 @@ private:
    string messageScfRmsDensity;
    string messageScfDampingThresh;
    string messageScfDampingWeight;
+   string messageScfDiisNumErrorVect;
+   string messageScfDiisStartError;
+   string messageScfDiisEndError;
    string stringSpace;
    string stringCommentOut;
    string stringTheory;
@@ -59,6 +62,9 @@ private:
    string stringScfRmsDensity;
    string stringScfDampingThresh;
    string stringScfDampingWeight;
+   string stringScfDiisNumErrorVect;
+   string stringScfDiisStartError;
+   string stringScfDiisEndError;
    string stringInertiaTensor;
    string stringInertiaTensorEnd;
    string stringInertiaTensorOrigin;
@@ -100,6 +106,9 @@ InputParser::InputParser(){
    this->messageScfRmsDensity = "\t\tRMS density: ";
    this->messageScfDampingThresh = "\t\tDamping threshold: ";
    this->messageScfDampingWeight = "\t\tDamping weight: ";
+   this->messageScfDiisNumErrorVect = "\t\tDIIS number of error vectors: ";
+   this->messageScfDiisStartError = "\t\tDIIS starting error: ";
+   this->messageScfDiisEndError = "\t\tDIIS ending error: ";
    this->stringSpace = " ";
    this->stringCommentOut = "//";
    this->stringTheoryCNDO2 = "cndo/2";
@@ -119,6 +128,9 @@ InputParser::InputParser(){
    this->stringScfRmsDensity = "rms_density";
    this->stringScfDampingThresh = "damping_thresh";
    this->stringScfDampingWeight = "damping_weight";
+   this->stringScfDiisNumErrorVect = "diis_num_error_vect";
+   this->stringScfDiisStartError = "diis_start_error";
+   this->stringScfDiisEndError = "diis_end_error";
    this->stringInertiaTensor = "inertia";
    this->stringInertiaTensorEnd = "inertia_end";
    this->stringInertiaTensorOrigin = "origin";
@@ -261,6 +273,21 @@ void InputParser::Parse(Molecule* molecule){
                Parameters::GetInstance()->SetDampingWeightSCF(atof(inputTerms[j+1].c_str()));
                j++;
             }
+            // DIIS number of stored error vectors
+            if(inputTerms[j].compare(this->stringScfDiisNumErrorVect) == 0){
+               Parameters::GetInstance()->SetDiisNumErrorVectSCF(atoi(inputTerms[j+1].c_str()));
+               j++;
+            }
+            // DIIS starting error
+            if(inputTerms[j].compare(this->stringScfDiisStartError) == 0){
+               Parameters::GetInstance()->SetDiisStartErrorSCF(atof(inputTerms[j+1].c_str()));
+               j++;
+            }
+            // DIIS ending error
+            if(inputTerms[j].compare(this->stringScfDiisEndError) == 0){
+               Parameters::GetInstance()->SetDiisEndErrorSCF(atof(inputTerms[j+1].c_str()));
+               j++;
+            }
             j++;   
          }
          i = j;
@@ -456,6 +483,9 @@ void InputParser::OutputScfConditions(){
    printf("%s%e\n",this->messageScfRmsDensity.c_str(),Parameters::GetInstance()->GetThresholdSCF());
    printf("%s%e\n",this->messageScfDampingThresh.c_str(),Parameters::GetInstance()->GetDampingThreshSCF());
    printf("%s%e\n",this->messageScfDampingWeight.c_str(),Parameters::GetInstance()->GetDampingWeightSCF());
+   printf("%s%d\n",this->messageScfDiisNumErrorVect.c_str(),Parameters::GetInstance()->GetDiisNumErrorVectSCF());
+   printf("%s%e\n",this->messageScfDiisStartError.c_str(),Parameters::GetInstance()->GetDiisStartErrorSCF());
+   printf("%s%e\n",this->messageScfDiisEndError.c_str(),Parameters::GetInstance()->GetDiisEndErrorSCF());
    cout << "\n";
 
 }

--- a/src/base/Parameters.h
+++ b/src/base/Parameters.h
@@ -23,6 +23,12 @@ public:
    void SetDampingThreshSCF(double dampingThreshSCF);
    double GetDampingWeightSCF();
    void SetDampingWeightSCF(double dampingWeightSCF);
+   int GetDiisNumErrorVectSCF();
+   void SetDiisNumErrorVectSCF(int diisNumErrorVectSCF);
+   double GetDiisStartErrorSCF();
+   void SetDiisStartErrorSCF(double diisStartErrorSCF);
+   double GetDiisEndErrorSCF();
+   void SetDiisEndErrorSCF(double diisEndErrorSCF);
    double GetEV2AU();
    double GetAngstrom2AU();
    double GetKayser2AU();
@@ -62,6 +68,9 @@ private:
    int maxIterationsSCF;
    double dampingThreshSCF;
    double dampingWeightSCF;
+   int diisNumErrorVectSCF;
+   double diisStartErrorSCF;
+   double diisEndErrorSCF;
    void SetDefaultValues();
    double eV2AU;
    double angstrom2AU;
@@ -123,6 +132,9 @@ void Parameters::SetDefaultValues(){
    this->maxIterationsSCF = 100;
    this->dampingThreshSCF = 1.0;
    this->dampingWeightSCF = 0.8;
+   this->diisNumErrorVectSCF = 5;
+   this->diisStartErrorSCF = pow(10.0, -2.0);
+   this->diisEndErrorSCF = pow(10.0, -8.0);
    this->currentTheory = CNDO2;
    this->gMolin2AU = pow(10.0,5.0)/(6.0221415*9.1095);
    this->degree2Radian = M_PI / 180.0;
@@ -175,6 +187,30 @@ void Parameters::SetDampingWeightSCF(double dampingWeightSCF){
    this->dampingWeightSCF = dampingWeightSCF;
 }
 
+int Parameters::GetDiisNumErrorVectSCF(){
+   return this->diisNumErrorVectSCF;
+}
+
+void Parameters::SetDiisNumErrorVectSCF(int diisNumErrorVectSCF){
+   this->diisNumErrorVectSCF = diisNumErrorVectSCF;
+}
+
+double Parameters::GetDiisStartErrorSCF(){
+   return this->diisStartErrorSCF;
+}
+
+void Parameters::SetDiisStartErrorSCF(double diisStartErrorSCF){
+   this->diisStartErrorSCF = diisStartErrorSCF;
+}
+
+double Parameters::GetDiisEndErrorSCF(){
+   return this->diisEndErrorSCF;
+}
+
+void Parameters::SetDiisEndErrorSCF(double diisEndErrorSCF){
+   this->diisEndErrorSCF = diisEndErrorSCF;
+}
+
 double Parameters::GetEV2AU(){
    return this->eV2AU;
 }

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -102,6 +102,20 @@ private:
    double GetAuxiliaryA(int k, double rho);
    double GetAuxiliaryB(int k, double rho);
    double GetAuxiliaryD(int la, int lb, int m);
+   void DoesDamp(double rmsDensity, 
+                 double** orbitalElectronPopulation, 
+                 double** oldOrbitalElectronPopulation, 
+                 Molecule* molecule);
+   void DoesDIIS(double** orbitalElectronPopulation,
+                 double** oldOrbitalElectronPopulation,
+                 double* atomicElectronPopulation,
+                 double*** diisStoredDensityMatrix,
+                 double*** diisStoredErrorVect,
+                 double** diisErrorProducts,
+                 double* diisErrorCoefficients,
+                 int diisNumErrorVect,
+                 Molecule* molecule, 
+                 int step);
    void OutputResults(double** fockMatrix, double* energiesMO, double* atomicElectronPopulation, Molecule* molecule);
    void CheckEnableAtomType(Molecule* molecule);
    void CheckNumberValenceElectrons(Molecule* molecule);
@@ -111,6 +125,11 @@ private:
                                 double** hMatrix, 
                                 double** dammyOrbitalElectronPopulation, 
                                 double*  dammyAtomicElectronPopulation );
+   void FreeSCFTemporaryMatrices(double** oldOrbitalElectronPopulation,
+                                 double*** diisStoredDensityMatrix,
+                                 double*** diisStoredErrorVect,
+                                 double** diisErrorProducts,
+                                 double* diisErrorCoefficients);
 
 };
 
@@ -289,8 +308,28 @@ void Cndo2::DoesSCF(){
       throw MolDSException(ss.str());
    }
 
+   // diis parameters
+   int diisNumErrorVect = Parameters::GetInstance()->GetDiisNumErrorVectSCF();
+
+
+   // malloc temporary matrices for scf
    double** oldOrbitalElectronPopulation = MallocerFreer::GetInstance()->MallocDoubleMatrix2d
             (this->molecule->GetTotalNumberAOs(), this->molecule->GetTotalNumberAOs());
+
+   // malloc temporary matrices for diis
+   double*** diisStoredDensityMatrix = NULL;
+   double*** diisStoredErrorVect = NULL;
+   double** diisErrorProducts = NULL;
+   double* diisErrorCoefficients = NULL;
+   if(0<diisNumErrorVect){
+      diisStoredDensityMatrix = MallocerFreer::GetInstance()->MallocDoubleMatrix3d
+                (diisNumErrorVect, this->molecule->GetTotalNumberAOs(), this->molecule->GetTotalNumberAOs());
+      diisStoredErrorVect = MallocerFreer::GetInstance()->MallocDoubleMatrix3d
+                (diisNumErrorVect, this->molecule->GetTotalNumberAOs(), this->molecule->GetTotalNumberAOs());
+      diisErrorProducts = MallocerFreer::GetInstance()->MallocDoubleMatrix2d(diisNumErrorVect+1, diisNumErrorVect+1);
+      diisErrorCoefficients = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(diisNumErrorVect+1);
+   }
+
    try{
       // calculate electron integral
       this->CalcGammaAB(this->gammaAB, this->molecule);
@@ -298,8 +337,6 @@ void Cndo2::DoesSCF(){
 
       // SCF
       double rmsDensity;
-      double dampingWeight = Parameters::GetInstance()->GetDampingWeightSCF();
-      double dampingThresh = Parameters::GetInstance()->GetDampingThreshSCF();
       int maxIterationsSCF = Parameters::GetInstance()->GetMaxIterationsSCF();
       bool isGuess=true;
       for(int i=0; i<maxIterationsSCF; i++){
@@ -312,7 +349,7 @@ void Cndo2::DoesSCF(){
                               this->orbitalElectronPopulation, 
                               this->atomicElectronPopulation,
                               isGuess);
-   
+
          // diagonalization
          bool calcEigenVectors = true;
          MolDS_mkl_wrapper::LapackWrapper::GetInstance()->Dsyevd(this->fockMatrix, 
@@ -325,17 +362,6 @@ void Cndo2::DoesSCF(){
                                              this->molecule, 
                                              this->fockMatrix);
 
-         // damping
-         if(dampingThresh < rmsDensity){
-            for(int j=0; j<this->molecule->GetTotalNumberAOs(); j++){
-               for(int k=0; k<this->molecule->GetTotalNumberAOs(); k++){
-                  this->orbitalElectronPopulation[j][k] *= (1.0-dampingWeight);
-                  this->orbitalElectronPopulation[j][k] += dampingWeight*oldOrbitalElectronPopulation[j][k];
-                                       
-               }
-            }
-         }
-
          // calc. electron population in each atom.
          this->CalcAtomicElectronPopulation(this->atomicElectronPopulation, 
                                             this->orbitalElectronPopulation, 
@@ -347,13 +373,31 @@ void Cndo2::DoesSCF(){
                                               this->molecule->GetTotalNumberAOs(), &rmsDensity, i)){
 
             cout << this->messageSCFMetConvergence;
-            this->OutputResults(this->fockMatrix, 
-                            this->energiesMO, 
-                            this->atomicElectronPopulation, 
-                            this->molecule);
-
+            this->OutputResults(this->fockMatrix, this->energiesMO, this->atomicElectronPopulation, this->molecule);
             break;
          }
+         else{
+           
+            // damping
+            this->DoesDamp(rmsDensity, this->orbitalElectronPopulation, oldOrbitalElectronPopulation, this->molecule);
+           
+            // diis 
+            this->DoesDIIS(this->orbitalElectronPopulation,
+                           oldOrbitalElectronPopulation,
+                           this->atomicElectronPopulation,
+                           diisStoredDensityMatrix,
+                           diisStoredErrorVect,
+                           diisErrorProducts,
+                           diisErrorCoefficients,
+                           diisNumErrorVect,
+                           this->molecule,
+                           i);
+
+            this->UpdateOldOrbitalElectronPopulation(oldOrbitalElectronPopulation, 
+                                                     this->orbitalElectronPopulation, 
+                                                     this->molecule->GetTotalNumberAOs());
+
+         }
 
          // SCF fails
          if(i==maxIterationsSCF-1){
@@ -361,27 +405,169 @@ void Cndo2::DoesSCF(){
             ss << this->errorMessageSCFNotConverged << maxIterationsSCF << "\n";
             throw MolDSException(ss.str());
          }
-
       }
    }
    catch(MolDSException ex){
-      if(oldOrbitalElectronPopulation != NULL){
-         MallocerFreer::GetInstance()->FreeDoubleMatrix2d
-         (oldOrbitalElectronPopulation, this->molecule->GetTotalNumberAOs());
-         oldOrbitalElectronPopulation = NULL;
-         //cout << "oldOrbitalElectronPopulation deleted\n";
-      }
+      this->FreeSCFTemporaryMatrices(oldOrbitalElectronPopulation,
+                                     diisStoredDensityMatrix,
+                                     diisStoredErrorVect,
+                                     diisErrorProducts,
+                                     diisErrorCoefficients);
+
       throw ex;
    }
-   
+   this->FreeSCFTemporaryMatrices(oldOrbitalElectronPopulation,
+                                  diisStoredDensityMatrix,
+                                  diisStoredErrorVect,
+                                  diisErrorProducts,
+                                  diisErrorCoefficients);
+
+
+   cout << this->messageDoneSCF;
+
+}
+
+void Cndo2::FreeSCFTemporaryMatrices(double** oldOrbitalElectronPopulation,
+                                     double*** diisStoredDensityMatrix,
+                                     double*** diisStoredErrorVect,
+                                     double** diisErrorProducts,
+                                     double* diisErrorCoefficients){
+
+   int diisNumErrorVect = Parameters::GetInstance()->GetDiisNumErrorVectSCF();
    if(oldOrbitalElectronPopulation != NULL){
       MallocerFreer::GetInstance()->FreeDoubleMatrix2d
       (oldOrbitalElectronPopulation, this->molecule->GetTotalNumberAOs());
       oldOrbitalElectronPopulation = NULL;
       //cout << "oldOrbitalElectronPopulation deleted\n";
    }
+   if(diisStoredDensityMatrix != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix3d
+      (diisStoredDensityMatrix, diisNumErrorVect, this->molecule->GetTotalNumberAOs());
+      diisStoredDensityMatrix = NULL;
+      //cout << "diisStoredDensityMatrix deleted\n";
+   } 
+   if(diisStoredErrorVect != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix3d
+      (diisStoredErrorVect, diisNumErrorVect, this->molecule->GetTotalNumberAOs());
+      diisStoredErrorVect = NULL;
+      //cout << "diisStoredErrorVect deleted\n";
+   } 
+   if(diisErrorProducts != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix2d
+      (diisErrorProducts, diisNumErrorVect+1);
+      diisErrorProducts = NULL;
+      //cout << "diisErrorProducts deleted\n";
+   } 
+   if(diisErrorCoefficients != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d
+      (diisErrorCoefficients);
+      diisErrorCoefficients = NULL;
+      //cout << "diisErrorCoefficients deleted\n";
+   } 
 
-   cout << this->messageDoneSCF;
+}
+
+/***
+ *
+ *  see ref. [P_1980] for diis methods.
+ *
+ */
+void Cndo2::DoesDIIS(double** orbitalElectronPopulation,
+                     double** oldOrbitalElectronPopulation,
+                     double* atomicElectronPopulation,
+                     double*** diisStoredDensityMatrix,
+                     double*** diisStoredErrorVect,
+                     double** diisErrorProducts,
+                     double* diisErrorCoefficients,
+                     int diisNumErrorVect,
+                     Molecule* molecule,
+                     int step){
+
+   int totalNumberAOs = molecule->GetTotalNumberAOs();
+   double diisStartError = Parameters::GetInstance()->GetDiisStartErrorSCF();
+   double diisEndError = Parameters::GetInstance()->GetDiisEndErrorSCF();
+
+   if( 0 < diisNumErrorVect){
+      for(int m=0; m<diisNumErrorVect-1; m++){
+         for(int j=0; j<totalNumberAOs; j++){
+            for(int k=0; k<totalNumberAOs; k++){
+               diisStoredDensityMatrix[m][j][k] = diisStoredDensityMatrix[m+1][j][k];
+               diisStoredErrorVect[m][j][k] = diisStoredErrorVect[m+1][j][k];
+            }
+         }
+      }
+
+      for(int j=0; j<totalNumberAOs; j++){
+         for(int k=0; k<totalNumberAOs; k++){
+            diisStoredDensityMatrix[diisNumErrorVect-1][j][k] = orbitalElectronPopulation[j][k];
+            diisStoredErrorVect[diisNumErrorVect-1][j][k] = orbitalElectronPopulation[j][k] 
+                                                           -oldOrbitalElectronPopulation[j][k];
+                     
+         }
+      }
+
+      for(int mi=0; mi<diisNumErrorVect-1; mi++){
+         for(int mj=0; mj<diisNumErrorVect-1; mj++){
+            diisErrorProducts[mi][mj] = diisErrorProducts[mi+1][mj+1];
+         }
+      }
+               
+      double tempErrorProduct=0.0;
+      for(int mi=0; mi<diisNumErrorVect; mi++){
+         tempErrorProduct = 0.0;
+         for(int j=0; j<totalNumberAOs; j++){
+            for(int k=0; k<totalNumberAOs; k++){
+               tempErrorProduct += diisStoredErrorVect[mi][j][k]*diisStoredErrorVect[diisNumErrorVect-1][j][k];
+            }
+         }
+         diisErrorProducts[mi][diisNumErrorVect-1] = tempErrorProduct;
+         diisErrorProducts[diisNumErrorVect-1][mi] = tempErrorProduct;
+         diisErrorProducts[mi][diisNumErrorVect] = -1.0;
+         diisErrorProducts[diisNumErrorVect][mi] = -1.0;
+         diisErrorCoefficients[mi] = 0.0;
+      }
+      diisErrorProducts[diisNumErrorVect][diisNumErrorVect] = 0.0;
+      diisErrorCoefficients[diisNumErrorVect] = -1.0;
+
+      double eMax = 0;
+      for(int j=0; j<totalNumberAOs; j++){
+         for(int k=0; k<totalNumberAOs; k++){
+            eMax = max(eMax, fabs(diisStoredErrorVect[diisNumErrorVect-1][j][k]));
+         }
+      }
+
+      if(diisNumErrorVect <= step && diisEndError<eMax && eMax<diisStartError){
+         MolDS_mkl_wrapper::LapackWrapper::GetInstance()->Dsysv(diisErrorProducts, 
+                                                                diisErrorCoefficients, 
+                                                                diisNumErrorVect+1);
+         for(int j=0; j<totalNumberAOs; j++){
+            for(int k=0; k<totalNumberAOs; k++){
+               orbitalElectronPopulation[j][k] = 0.0;
+               for(int m=0; m<diisNumErrorVect; m++){
+                  orbitalElectronPopulation[j][k] += diisErrorCoefficients[m]*diisStoredDensityMatrix[m][j][k];
+               }
+            }
+         }
+
+         // calc. electron population in each atom.
+         this->CalcAtomicElectronPopulation(atomicElectronPopulation, orbitalElectronPopulation, molecule);
+      }
+   }
+}
+
+void Cndo2::DoesDamp(double rmsDensity, double** orbitalElectronPopulation, 
+                     double** oldOrbitalElectronPopulation, Molecule* molecule){
+
+   double dampingThresh = Parameters::GetInstance()->GetDampingThreshSCF();
+   double dampingWeight = Parameters::GetInstance()->GetDampingWeightSCF();
+   if(0.0 < dampingWeight && dampingThresh < rmsDensity){
+      for(int j=0; j<molecule->GetTotalNumberAOs(); j++){
+         for(int k=0; k<molecule->GetTotalNumberAOs(); k++){
+            orbitalElectronPopulation[j][k] *= (1.0 - dampingWeight);
+            orbitalElectronPopulation[j][k] += dampingWeight*oldOrbitalElectronPopulation[j][k];
+         }
+      }
+   } 
 
 }
 
@@ -596,11 +782,6 @@ bool Cndo2::SatisfyConvergenceCriterion(double** oldOrbitalElectronPopulation,
    if(*rmsDensity < Parameters::GetInstance()->GetThresholdSCF()){
       satisfy = true;
    }
-   else{
-      this->UpdateOldOrbitalElectronPopulation(oldOrbitalElectronPopulation, 
-                                               orbitalElectronPopulation, 
-                                               numberAOs);
-   }
 
    return satisfy; 
 }

--- a/src/input.in
+++ b/src/input.in
@@ -1,6 +1,6 @@
 THEORY
-   //cndo/2
-   indo
+   cndo/2
+   //indo
    //zindo/s
    //none 
    //principal_axes
@@ -9,15 +9,18 @@ THEORY
 THEORY_END
 
 SCF
-   max_iter 200
-   rms_density 0.00000001
+   max_iter 50
+   rms_density 0.000001
    damping_thresh 1.0 
-   damping_weight 0.7
+   damping_weight 0.0
+   diis_num_error_vect 5
+   diis_start_error 0.002
+   diis_end_error 0.00000002
 SCF_END
 
 CIS
-   activeOcc 100
-   activeVir 100
+   activeOcc 2 
+   activeVir 2
    nstates 1000
 CIS_END
 

--- a/src/mkl_wrapper/LapackWrapper.h
+++ b/src/mkl_wrapper/LapackWrapper.h
@@ -18,24 +18,32 @@ public:
    static LapackWrapper* GetInstance();
    static void DeleteInstance();
    int Dsyevd(double** matrix, double* eigenValues, int size, bool calcEigenVectors);
+   int Dsysv(double** matrix, double* b, int size);
 private:
    LapackWrapper();
    LapackWrapper(LapackWrapper&);
    void operator = (LapackWrapper&);
    ~LapackWrapper();
    static LapackWrapper* lapackWrapper;
-
+   bool calculatedDsysvWorkSize;
+   int dsysvWorkSize;
    string errorMessageDsyevdInfo;
    string errorMessageDsyevdSize;
+   string errorMessageDsysvInfo;
+   string errorMessageDsysvSize;
 };
 LapackWrapper* LapackWrapper::lapackWrapper = NULL;
 
 LapackWrapper::LapackWrapper(){
    this->errorMessageDsyevdInfo = "Error in mkl_wrapper::LapackWrapper::Dsyevd: info != 0\n";
    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";
 }
 
 LapackWrapper::~LapackWrapper(){
+   this->calculatedDsysvWorkSize = false;
+   this->dsysvWorkSize = 64;
 }
 
 LapackWrapper* LapackWrapper::GetInstance(){
@@ -154,6 +162,72 @@ int LapackWrapper::Dsyevd(double** matrix, double* eigenValues, int size, bool c
    return info;
 }
 
+/***
+ *
+ * Slove matrix*X=b, then we get X by this method.
+ * The X is stored in b.
+ *
+ */
+int LapackWrapper::Dsysv(double** matrix, double* b, int size){
+   int info = 0;
+   int lwork;
+   char uplo = 'U';
+   int lda = size;
+   int ldb = size;
+   int nrhs = 1;
+   double* convertedMatrix;
+   double* work;
+   int* ipiv;
+
+   if(size < 1 ){
+      stringstream ss;
+      ss << errorMessageDsysvSize;
+      throw MolDSException(ss.str());
+   }
+
+   // malloc
+   ipiv = MallocerFreer::GetInstance()->MallocIntMatrix1d(2*size);
+   convertedMatrix = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(size*size);
+
+   for(int i = 0; i < size; i++){
+      for(int j = i; j < size; j++){
+         convertedMatrix[i+j*size] = matrix[i][j];
+      }
+   }
+
+   // calc. lwork
+   if(!this->calculatedDsysvWorkSize){
+      lwork = -1;
+      double tempWork[3]={0.0, 0.0, 0.0};
+      dsysv(&uplo, &size, &nrhs, convertedMatrix, &lda, ipiv, b, &ldb, tempWork, &lwork, &info);
+      this->calculatedDsysvWorkSize = true;
+      this->dsysvWorkSize = tempWork[0];
+   }
+
+   info = 0;
+   lwork = this->dsysvWorkSize;
+   work = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(lwork);
+
+   // call Lapack
+   dsysv(&uplo, &size, &nrhs, convertedMatrix, &lda, ipiv, b, &ldb, work, &lwork, &info);
+
+   // free
+   MallocerFreer::GetInstance()->FreeDoubleMatrix1d(convertedMatrix);
+   convertedMatrix = NULL;
+   MallocerFreer::GetInstance()->FreeIntMatrix1d(ipiv);
+   ipiv = NULL;
+   MallocerFreer::GetInstance()->FreeDoubleMatrix1d(work);
+   work = NULL;
+  
+   if(info != 0){
+      cout << "info=" << info << endl;
+      stringstream ss;
+      ss << errorMessageDsysvInfo;
+      throw MolDSException(ss.str());
+   }
+   return info;
+}
+
 
 
 }