提交: 5998581cb936ac699fdf1f9e987b8f72cdea22c6 - MolDS (git)

Commit MetaInfo

• 修订版: 5998581cb936ac699fdf1f9e987b8f72cdea22c6 (tree)
• 时间: 2011-01-05 18:59:57
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Inertia tensor is calculated.

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

更改概述

• modified: doc/README.txt
• modified: src/MolDS.cpp
• modified: src/base/Enums.h
• modified: src/base/InputParser.h
• modified: src/base/Molecule.h
• modified: src/base/Parameters.h
• modified: src/base/atoms/Atom.h
• modified: src/base/atoms/Catom.h
• modified: src/base/atoms/Hatom.h
• modified: src/base/atoms/Liatom.h
• modified: src/base/atoms/Satom.h
• modified: src/input.in

差异

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -1,5 +1,30 @@
-for 32 bit
+
+Compile: 
+   for 32 bit
    $icc MolDS.cpp -lmkl_intel -lmkl_intel_thread -lmkl_core -liomp5 -lpthread
 
-for 64 bit
+   for 64 bit
    $icc MolDS.cpp -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread
+
+
+Carring Out:
+   $./a.out < input.in
+
+Principal Axes (Diagonalizing the inertia tensor):
+   Write "principal-axes" in theory-directive.
+
+   E.g. 
+         THEORY
+            principal-axes
+         THEORY_END
+   
+   -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.
+
+    E.g.
+         INERTIA
+            origin 1.2 2.3 3.4
+         INERTIA_END
+

--- a/src/MolDS.cpp
+++ b/src/MolDS.cpp
@@ -108,6 +108,17 @@ int main(){
          delete zindoS;
       }
 
+      // Diagonalize Inertia Tensor
+      else if(Parameters::GetInstance()->GetCurrentTheory() == PrincipalAxes && runingNormally){
+         try{
+            molecule->CalcPrincipalAxes();
+         }
+         catch(MolDSException ex){
+            cout << ex.what() << endl;
+            runingNormally = false;
+         }
+
+      }
 
       //Free 
       LapackWrapper::DeleteInstance(); 

--- a/src/base/Enums.h
+++ b/src/base/Enums.h
@@ -11,6 +11,7 @@ RENUMSTR_BEGIN( TheoryType, TheoryTypeStr )
    RENUMSTR( CNDO2,  "CNDO/2" )
    RENUMSTR( INDO,   "INDO" )
    RENUMSTR( ZINDOS, "ZINDO/S" )
+   RENUMSTR( PrincipalAxes, "Principal Axes" )
    RENUMSTR( NONE, "NONE" )
    RENUMSTR( TheoryType_end,  "TheoryType_end" )
 RENUMSTR_END()

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -43,11 +43,15 @@ private:
    string stringTheoryCNDO2;
    string stringTheoryINDO;
    string stringTheoryZINDOS;
+   string stringTheoryPrincipalAxes;
    string stringTheoryNONE;
    string stringScf;
    string stringScfEnd;
    string stringScfMaxIter;
    string stringScfRmsDensity;
+   string stringInertiaTensor;
+   string stringInertiaTensorEnd;
+   string stringInertiaTensorOrigin;
    void CalcMolecularBasics(Molecule* molecule);
    void OutputMolecularBasics(Molecule* molecule);
    void OutputScfConditions();
@@ -86,6 +90,10 @@ InputParser::InputParser(){
    this->stringTheoryCNDO2 = "cndo/2";
    this->stringTheoryINDO = "indo";
    this->stringTheoryZINDOS = "zindo/s";
+   this->stringTheoryPrincipalAxes = "principal-axes";
+   this->stringInertiaTensor = "inertia";
+   this->stringInertiaTensorEnd = "inertia_end";
+   this->stringInertiaTensorOrigin = "origin";
    this->stringTheoryNONE = "none";
 }
 
@@ -159,9 +167,9 @@ void InputParser::Parse(Molecule* molecule){
       if(inputTerms[i].compare(this->stringGeometry) == 0){
          int j=i+1;
          while(inputTerms[j].compare(this->stringGeometryEnd) != 0){
-            double x = atof(inputTerms[j+1].c_str());
-            double y = atof(inputTerms[j+2].c_str());
-            double z = atof(inputTerms[j+3].c_str());
+            double x = atof(inputTerms[j+1].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
+            double y = atof(inputTerms[j+2].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
+            double z = atof(inputTerms[j+3].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
             if(inputTerms[j] == "h"){
                molecule->GetAtomVect()->push_back(new Hatom(x, y, z));
             }
@@ -198,6 +206,23 @@ void InputParser::Parse(Molecule* molecule){
          i = j;
       }
       
+      // inertia tensor condition
+      if(inputTerms[i].compare(this->stringInertiaTensor) == 0){
+         int j=i+1;
+         while(inputTerms[j].compare(this->stringInertiaTensorEnd) != 0){
+            // origin
+            if(inputTerms[j].compare(this->stringInertiaTensorOrigin) == 0){
+               double x = atof(inputTerms[j+1].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
+               double y = atof(inputTerms[j+2].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
+               double z = atof(inputTerms[j+3].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
+               molecule->SetInertiaTensorOrigin(x, y, z);
+               j+=3;
+            }
+            j++;   
+         }
+         i = j;
+      }
+      
       // theory
       if(inputTerms[i].compare(this->stringTheory) == 0){
          int j=i+1;
@@ -218,10 +243,16 @@ void InputParser::Parse(Molecule* molecule){
                Parameters::GetInstance()->SetCurrentTheory(ZINDOS);
             }
 
+            // Princepal axes
+            else if(inputTerms[j].compare(this->stringTheoryPrincipalAxes) == 0){
+               Parameters::GetInstance()->SetCurrentTheory(PrincipalAxes);
+            }
+
             // NONE 
             else if(inputTerms[j].compare(this->stringTheoryNONE) == 0){
                Parameters::GetInstance()->SetCurrentTheory(NONE);
             }
+
             j++;
          }
          i = j;
@@ -300,6 +331,7 @@ bool InputParser::IsCommentOut(string tempStr){
 
 }
 
+
 }
 #endif
 

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -16,9 +16,12 @@ class Molecule{
 private:
    vector<Atom*>* atomVect;
    double* COMXyz;
+   double* inertiaTensorOrigin;
    bool wasCalculatedCOMXyz;
    int totalNumberAOs;
    int totalNumberValenceElectrons;
+   void CalcInertiaTensor(double** inertiaTensor);
+   void FreeInertiaTensorMoments(double** inertiaTensor, double* inertiaMoments);
    string messageTotalNumberAOs;
    string messageTotalNumberAtoms;
    string messageTotalNumberValenceElectrons;
@@ -41,11 +44,14 @@ public:
    void OutputCOMXyz();
    void OutputTotalNumberAtomsAOsValenceelectrons();
    void OutputConfiguration();
+   void SetInertiaTensorOrigin(double x, double y, double z);
+   void CalcPrincipalAxes();
 };
 
 Molecule::Molecule(){
    this->atomVect = new vector<Atom*>;
    this->COMXyz = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
+   this->inertiaTensorOrigin = NULL;
    this->wasCalculatedCOMXyz = false;
    this->messageTotalNumberAOs = "\tTotal number of valence AOs: ";
    this->messageTotalNumberAtoms = "\tTotal number of atoms: ";
@@ -73,6 +79,11 @@ Molecule::~Molecule(){
       this->COMXyz = NULL;
       //cout << "COMXyz deleted\n";
    }
+   if(this->inertiaTensorOrigin != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(this->inertiaTensorOrigin);
+      this->inertiaTensorOrigin = NULL;
+      //cout << "inertiaTensorOrigin deleted\n";
+   }
 }
 
 vector<Atom*>* Molecule::GetAtomVect(){
@@ -175,6 +186,89 @@ void Molecule::OutputTotalNumberAtomsAOsValenceelectrons(){
    cout << this->messageTotalNumberValenceElectrons << this->totalNumberValenceElectrons << "\n\n";
 }
 
+void Molecule::SetInertiaTensorOrigin(double x, double y, double z){
+   if(this->inertiaTensorOrigin == NULL){
+      this->inertiaTensorOrigin = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
+   }
+
+   this->inertiaTensorOrigin[0] = x;
+   this->inertiaTensorOrigin[1] = y;
+   this->inertiaTensorOrigin[2] = z;
+
+}
+
+void Molecule::CalcPrincipalAxes(){
+
+   if(this->inertiaTensorOrigin == NULL){
+      this->SetInertiaTensorOrigin(this->COMXyz[0], this->COMXyz[1], this->COMXyz[2]);
+   }
+
+   double** inertiaTensor = MallocerFreer::GetInstance()->MallocDoubleMatrix2d(3, 3);
+   double*  inertiaMoments = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
+
+   try{
+      this->CalcInertiaTensor(inertiaTensor);
+      
+      // ToDo: diagonalization!!!!!!
+      
+   }
+   catch(MolDSException ex){
+      this->FreeInertiaTensorMoments(inertiaTensor, inertiaMoments);
+      throw ex;
+   }
+
+   this->FreeInertiaTensorMoments(inertiaTensor, inertiaMoments);
+   
+}
+
+void Molecule::CalcInertiaTensor(double** inertiaTensor){
+
+   Atom* atom;
+   double x;
+   double y;
+   double z;
+   double atomicMass;
+   for(int a=0; a<this->atomVect->size(); a++){
+      atom = (*this->atomVect)[a];
+      atomicMass = atom->GetAtomicMass();
+      x = atom->GetXyz()[0] - this->inertiaTensorOrigin[0];
+      y = atom->GetXyz()[1] - this->inertiaTensorOrigin[1];
+      z = atom->GetXyz()[2] - this->inertiaTensorOrigin[2];
+
+      inertiaTensor[0][0] += atomicMass*(y*y + z*z);
+      inertiaTensor[0][1] -= atomicMass*x*y;
+      inertiaTensor[0][2] -= atomicMass*x*z;
+
+      inertiaTensor[1][0] -= atomicMass*y*x;
+      inertiaTensor[1][1] += atomicMass*(x*x + z*z);
+      inertiaTensor[1][2] -= atomicMass*y*z;
+
+      inertiaTensor[2][0] -= atomicMass*z*x;
+      inertiaTensor[2][1] -= atomicMass*z*y;
+      inertiaTensor[2][2] += atomicMass*(x*x + y*y);
+
+   }
+   
+}
+
+void Molecule::FreeInertiaTensorMoments(double** inertiaTensor, double* inertiaMoments){
+
+   if(inertiaTensor != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix2d(inertiaTensor, 3);
+      inertiaTensor = NULL;
+      //cout << "inertiaTensor deleted\n";
+   }
+
+   if(inertiaMoments != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(inertiaMoments);
+      inertiaMoments = NULL;
+      //cout << "inertiaMoments deleted\n";
+   }
+
+}
+
+
+
 }
 #endif
 

--- a/src/base/Parameters.h
+++ b/src/base/Parameters.h
@@ -24,6 +24,7 @@ private:
    double eV2AU;
    double angstrom2AU;
    double kayser2AU;
+   double gMolin2AU;
    double bondingAdjustParameterK; //see (3.79) in J. A. Pople book
    TheoryType currentTheory;
 
@@ -38,6 +39,7 @@ public:
    double GetEV2AU();
    double GetAngstrom2AU();
    double GetKayser2AU();
+   double GetGMolin2AU();
    double GetBondingAdjustParameterK();
    TheoryType GetCurrentTheory();
    void SetCurrentTheory(TheoryType theory);
@@ -73,6 +75,7 @@ void Parameters::SetDefaultValues(){
    this->thresholdSCF = pow(10.0, -8.0);
    this->maxIterationsSCF = 100;
    this->currentTheory = CNDO2;
+   this->gMolin2AU = pow(10.0,5.0)/(6.0221415*9.1095);
 }
 
 double Parameters::GetThresholdSCF(){
@@ -103,6 +106,10 @@ double Parameters::GetKayser2AU(){
    return this->kayser2AU;
 }
 
+double Parameters::GetGMolin2AU(){
+   return this->gMolin2AU;
+}
+
 double Parameters::GetBondingAdjustParameterK(){
    return this->bondingAdjustParameterK;
 }

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -160,9 +160,9 @@ double* Atom::GetXyz(){
 }
 
 void Atom::SetXyz(double x, double y, double z){
-   xyz[0]= x * Parameters::GetInstance()->GetAngstrom2AU();
-   xyz[1]= y * Parameters::GetInstance()->GetAngstrom2AU();
-   xyz[2]= z * Parameters::GetInstance()->GetAngstrom2AU();
+   xyz[0]= x;
+   xyz[1]= y;
+   xyz[2]= z;
 }
 
 vector<OrbitalType> Atom::GetValence(){

--- a/src/base/atoms/Catom.h
+++ b/src/base/atoms/Catom.h
@@ -16,7 +16,7 @@ public:
 
 Catom::Catom(double x, double y, double z) : Atom(x, y, z){
    this->atomType = C;
-   this->atomicMass = 12.0107;
+   this->atomicMass = 12.0107*Parameters::GetInstance()->GetGMolin2AU();
    this->valence.push_back(s);
    this->valence.push_back(py);
    this->valence.push_back(pz);

--- a/src/base/atoms/Hatom.h
+++ b/src/base/atoms/Hatom.h
@@ -17,7 +17,7 @@ public:
 
 Hatom::Hatom(double x, double y, double z) : Atom(x, y, z){
    this->atomType = H;
-   this->atomicMass = 1.00794;
+   this->atomicMass = 1.00794*Parameters::GetInstance()->GetGMolin2AU();
    this->valence.push_back(s);
    this->bondingParameter = -9.0*Parameters::GetInstance()->GetEV2AU();
    this->bondingParameterSZindo = -9.0*Parameters::GetInstance()->GetEV2AU();

--- a/src/base/atoms/Liatom.h
+++ b/src/base/atoms/Liatom.h
@@ -17,7 +17,7 @@ public:
 
 Liatom::Liatom(double x, double y, double z) : Atom(x, y, z){
    this->atomType = Li;
-   this->atomicMass = 6.941;
+   this->atomicMass = 6.941*Parameters::GetInstance()->GetGMolin2AU();
    this->valence.push_back(s);
    this->valence.push_back(py);
    this->valence.push_back(pz);

--- a/src/base/atoms/Satom.h
+++ b/src/base/atoms/Satom.h
@@ -16,7 +16,7 @@ public:
 
 Satom::Satom(double x, double y, double z) : Atom(x, y, z){
    this->atomType = S;
-   this->atomicMass = 32.066;
+   this->atomicMass = 32.066*Parameters::GetInstance()->GetGMolin2AU();
    this->valence.push_back(s);
    this->valence.push_back(py);
    this->valence.push_back(pz);

--- a/src/input.in
+++ b/src/input.in
@@ -7,9 +7,15 @@ THEORY
    //cndo/2
    //indo
    //zindo/s
-   none 
+   //none 
+   principal-axes
 THEORY_END
 
+INERTIA
+   origin 1.0 2.0 3.0
+INERTIA_END
+
+
 //metane
 //GEOMETRY
 //   C  -0.37687006  0.95490165  0.00000000