提交: 1ec4a1fe9717c2b2a7dcc04f409287c5e5198c12 - MolDS (git)

Commit MetaInfo

• 修订版: 1ec4a1fe9717c2b2a7dcc04f409287c5e5198c12 (tree)
• 时间: 2011-10-19 22:14:12
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Heats of formation in MNDO is implemented. #26393

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

更改概述

• modified: src/cndo/Cndo2.h
• modified: src/input.in
• modified: src/mndo/Mndo.h

差异

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -46,6 +46,7 @@ protected:
    string messageUnitSec; 
    vector<AtomType> enableAtomTypes;
    double coreRepulsionEnergy;
+   virtual void CalcHFProperties();
    virtual void CalcCoreRepulsionEnergy();
    virtual double GetDiatomCoreRepulsionFirstDerivative(int indexAtomA, 
                                                         int indexAtomB, 
@@ -87,6 +88,10 @@ protected:
                                               Molecule* molecule, double** fockMatrix, double** gammaAB);
    virtual void CalcTwoElecTowCore(double****** twoElecTwoCore, Molecule* molecule);
    virtual void CalcForce(int electronicStateIndex);
+   virtual void OutputHFResults(double** fockMatrix, 
+                                double* energiesMO, 
+                                double* atomicElectronPopulation, 
+                                Molecule* molecule);
    TheoryType theory;
    Molecule* molecule;
    double** fockMatrix;
@@ -189,8 +194,6 @@ private:
                  int diisNumErrorVect,
                  Molecule* molecule, 
                  int step);
-   void OutputResults(double** fockMatrix, double* energiesMO, 
-                      double* atomicElectronPopulation, Molecule* molecule);
    void CheckEnableAtomType(Molecule* molecule);
    void CheckNumberValenceElectrons(Molecule* molecule);
    void FreeDiatomicOverlapAndRotatingMatrix(double*** diatomicOverlap, 
@@ -505,22 +508,11 @@ void Cndo2::DoesSCF(bool requiresGuess){
 
             // calc. some properties.
             // e.g. electronic energy, electron population in each atom, and core replsion.
-            this->CalcAtomicElectronPopulation(this->atomicElectronPopulation, 
-                                               this->orbitalElectronPopulation, 
-                                               this->molecule);
-
-            this->CalcCoreRepulsionEnergy();
-            this->CalcElecEnergy(&this->elecEnergy, 
-                                 this->molecule, 
-                                 this->energiesMO, 
-                                 this->fockMatrix, 
-                                 this->gammaAB,
-                                 this->coreRepulsionEnergy);
-
-            this->OutputResults(this->fockMatrix, 
-                                this->energiesMO, 
-                                this->atomicElectronPopulation, 
-                                this->molecule);
+            this->CalcHFProperties();
+            this->OutputHFResults(this->fockMatrix, 
+                                  this->energiesMO, 
+                                  this->atomicElectronPopulation, 
+                                  this->molecule);
             break;
          }
          else{
@@ -577,6 +569,19 @@ void Cndo2::DoesSCF(){
    this->DoesSCF(requiresGuess);
 }
 
+void Cndo2::CalcHFProperties(){
+   this->CalcAtomicElectronPopulation(this->atomicElectronPopulation, 
+                                      this->orbitalElectronPopulation, 
+                                      this->molecule);
+   this->CalcCoreRepulsionEnergy();
+   this->CalcElecEnergy(&this->elecEnergy, 
+                        this->molecule, 
+                        this->energiesMO, 
+                        this->fockMatrix, 
+                        this->gammaAB,
+                        this->coreRepulsionEnergy);
+}
+
 void Cndo2::DoesCIS(){
    stringstream ss;
    ss << this->errorMessageCISNotImplemented;
@@ -743,7 +748,7 @@ void Cndo2::DoesDamp(double rmsDensity, double** orbitalElectronPopulation,
 
 }
 
-void Cndo2::OutputResults(double** fockMatrix, double* energiesMO, double* atomicElectronPopulation, Molecule* molecule){
+void Cndo2::OutputHFResults(double** fockMatrix, double* energiesMO, double* atomicElectronPopulation, Molecule* molecule){
 
    // output MO energy
    cout << this->messageEnergiesMOs;
@@ -759,7 +764,7 @@ void Cndo2::OutputResults(double** fockMatrix, double* energiesMO, double* atomi
          mo, this->messageUnOcc.c_str(), energiesMO[mo], energiesMO[mo]/eV2AU);
       }
    }
-   cout << endl << endl;
+   cout << endl;
 
    // output total energy
    cout << this->messageElecEnergy;
@@ -781,7 +786,7 @@ void Cndo2::OutputResults(double** fockMatrix, double* energiesMO, double* atomi
       Atom* atom = (*molecule->GetAtomVect())[a];
       printf("\t\t%d\t%s\t%e\t%e\n",a,AtomTypeStr(atom->GetAtomType()),atom->GetCoreCharge(),atom->GetCoreCharge()-atomicElectronPopulation[a]);
    }
-   cout << endl << endl;
+   cout << endl;
 
 
 }

--- a/src/input.in
+++ b/src/input.in
@@ -20,15 +20,15 @@ SCF
    diis_end_error 0.00000002
 SCF_END
 
-//CIS
-//   davidson yes
-//   active_occ 100
-//   active_vir 100
-//   nstates 5
-//   max_iter 200
-//   max_dim 100
-//   norm_tol 0.000001
-//CIS_END
+CIS
+   davidson yes
+   active_occ 100
+   active_vir 100
+   nstates 5
+   max_iter 200
+   max_dim 100
+   norm_tol 0.000001
+CIS_END
 
 //MD
 //   total_steps 100
@@ -170,28 +170,28 @@ SCF_END
 //GEOMETRY_END
 
 // h2s-2
-//GEOMETRY
-//   S                  0.0000     0.0000     0.0000
-//   H                  1.2993     0.0000     0.0000
-//   H                 -0.1814     1.2865     0.0000
-//GEOMETRY_END
-
-// benzene
 GEOMETRY
-   C  -0.505391      0.552561      0.000000
-   C  0.889769      0.552561      0.000000
-   C  1.587307      1.760312      0.000000
-   C  0.889653      2.968821     -0.001199
-   C  -0.505172      2.968743     -0.001678
-   C  -1.202773      1.760537     -0.000682
-   H  -1.055150     -0.399756      0.000450
-   H  1.439277     -0.399952      0.001315
-   H  2.686987      1.760392      0.000634
-   H  1.439853      3.920964     -0.001258
-   H  -1.055294      3.921024     -0.002631
-   H  -2.302377      1.760720     -0.000862
+   S                  0.0000     0.0000     0.0000
+   H                  1.2993     0.0000     0.0000
+   H                 -0.1814     1.2865     0.0000
 GEOMETRY_END
 
+// benzene
+//GEOMETRY
+//   C  -0.505391      0.552561      0.000000
+//   C  0.889769      0.552561      0.000000
+//   C  1.587307      1.760312      0.000000
+//   C  0.889653      2.968821     -0.001199
+//   C  -0.505172      2.968743     -0.001678
+//   C  -1.202773      1.760537     -0.000682
+//   H  -1.055150     -0.399756      0.000450
+//   H  1.439277     -0.399952      0.001315
+//   H  2.686987      1.760392      0.000634
+//   H  1.439853      3.920964     -0.001258
+//   H  -1.055294      3.921024     -0.002631
+//   H  -2.302377      1.760720     -0.000862
+//GEOMETRY_END
+
 // acetylene
 //GEOMETRY
 // C                     0.24933   0.0   0.0 

--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -19,6 +19,11 @@ protected:
    virtual void SetMessages();
    virtual void SetEnableAtomTypes();
    virtual void CalcCoreRepulsionEnergy();
+   virtual void CalcHFProperties();
+   virtual void OutputHFResults(double** fockMatrix, 
+                                double* energiesMO, 
+                                double* atomicElectronPopulation, 
+                                Molecule* molecule);
    virtual double GetFockDiagElement(Atom* atomA, 
                                      int atomAIndex, 
                                      int mu, 
@@ -69,6 +74,10 @@ private:
    string errorMessageCalcTwoElecTwoCoreDiatomicNullMatrix;
    string errorMessageCalcTwoElecTwoCoreNullMatrix;
    string errorMessageCalcTwoElecTwoCoreDiatomicSameAtoms;
+   string messageHeatsFormation;
+   string messageHeatsFormationTitle;
+   double heatsFormation;
+   void CalcHeatsFormation(double* heatsFormation, Molecule* molecule);
    double GetElectronCoreAttraction(int atomAIndex, int atomBIndex, 
                                     int mu, int nu, double****** twoElecTwoCore);
    void CalcTwoElecTwoCoreDiatomic(double**** matrix, int atomAIndex, int atomBIndex);
@@ -88,6 +97,7 @@ Mndo::Mndo() : MolDS_zindo::ZindoS(){
    this->theory = MNDO;
    this->SetMessages();
    this->SetEnableAtomTypes();
+   this->heatsFormation = 0.0;
    //cout << "Mndo created\n";
 }
 
@@ -140,6 +150,8 @@ void Mndo::SetMessages(){
    this->messageSCFMetConvergence = "\n\n\n\t\tMNDO/S-SCF met convergence criterion(^^b\n\n\n";
    this->messageStartSCF = "**********  START: MNDO/S-SCF  **********\n";
    this->messageDoneSCF = "**********  DONE: MNDO/S-SCF  **********\n\n\n";
+   this->messageHeatsFormation = "\tHeats of formation:\n";
+   this->messageHeatsFormationTitle = "\t\t| [a.u.] | [Kcal/mol] | \n";
    this->messageStartCIS = "**********  START: MNDO/S-CIS  **********\n";
    this->messageDoneCIS = "**********  DONE: MNDO/S-CIS  **********\n\n\n";
    this->messageDavidsonConverge = "\n\n\t\tDavidson for MNDO-CIS met convergence criterion(^^b\n\n\n";
@@ -189,6 +201,33 @@ void Mndo::CalcCoreRepulsionEnergy(){
    this->coreRepulsionEnergy = energy;
 }
 
+void Mndo::CalcHeatsFormation(double* heatsFormation, Molecule* molecule){
+   *heatsFormation = this->GetElectronicEnergy();
+   for(int A=0; A<molecule->GetAtomVect()->size(); A++){
+      Atom* atom = (*molecule->GetAtomVect())[A];
+      *heatsFormation -= atom->GetMndoElecEnergyAtom();
+      *heatsFormation += atom->GetMndoHeatsFormAtom();
+   }
+}
+
+void Mndo::CalcHFProperties(){
+   MolDS_cndo::Cndo2::CalcHFProperties();
+   this->CalcHeatsFormation(&this->heatsFormation, this->molecule);
+}
+
+void Mndo::OutputHFResults(double** fockMatrix, double* energiesMO, 
+                           double* atomicElectronPopulation, Molecule* molecule){
+   MolDS_cndo::Cndo2::OutputHFResults(fockMatrix, 
+                                      energiesMO, 
+                                      atomicElectronPopulation, 
+                                      molecule);
+   // output heats of formation
+   cout << this->messageHeatsFormation;
+   cout << this->messageHeatsFormationTitle;
+   printf("\t\t%e\t%e\n\n",this->heatsFormation,
+                           this->heatsFormation/Parameters::GetInstance()->GetKcalMolin2AU());
+}
+
 double Mndo::GetFockDiagElement(Atom* atomA, int atomAIndex, int mu, 
                                  Molecule* molecule, double** gammaAB,
                                  double** orbitalElectronPopulation,