Commit: 8e1b5afb7d3da25153560595abb9b7e062808855 - MolDS (git)

Commit MetaInfo

• Revisão: 8e1b5afb7d3da25153560595abb9b7e062808855 (tree)
• Hora: 2011-01-07 20:14:24
• Autor: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Mensagem de Log

Logic to rotate molecule is added.

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

Mudança Sumário

• modified: src/MolDS.cpp
• modified: src/base/Enums.h
• modified: src/base/MathUtilities.h
• modified: src/base/Molecule.h
• modified: src/input.in

Diff

--- a/src/MolDS.cpp
+++ b/src/MolDS.cpp
@@ -13,6 +13,7 @@
 #undef INCLUDED_ENUMS
 #define RENUMSTR_BODY 1  
 #include"base/Enums.h"
+#include"base/MathUtilities.h"
 #include"base/EularAngle.h"
 #include"base/Molecule.h"
 #include"base/atoms/Atom.h"

--- a/src/base/Enums.h
+++ b/src/base/Enums.h
@@ -18,13 +18,18 @@ RENUMSTR_BEGIN( TheoryType, TheoryTypeStr )
    RENUMSTR( TheoryType_end,  "TheoryType_end" )
 RENUMSTR_END()
 
-
 RENUMSTR_BEGIN( RotatingType, RotatingTypeStr )
    RENUMSTR( Axis,  "Axis" )
    RENUMSTR( Eular,  "EularAngle" )
    RENUMSTR( RotatingType_end,  "RotatingType_end" )
 RENUMSTR_END()
 
+RENUMSTR_BEGIN( RotatedObjectType, RotatedObjectTypeStr )
+   RENUMSTR( System,  "System" )
+   RENUMSTR( Frame,  "Frame" )
+   RENUMSTR( RotatedObjectType_end,  "RotatedObjectType_end" )
+RENUMSTR_END()
+
 RENUMSTR_BEGIN( ShellType, ShellTypeStr )
    RENUMSTR( k,  "k" )
    RENUMSTR( l,  "l" )
@@ -32,6 +37,13 @@ RENUMSTR_BEGIN( ShellType, ShellTypeStr )
    RENUMSTR( ShellType_end,  "ShellType_end" )
 RENUMSTR_END()
 
+RENUMSTR_BEGIN( CartesianType, CartesianTypeStr )
+   RENUMSTR( XAxis,  "XAxis" )
+   RENUMSTR( YAxis,  "YAxis" )
+   RENUMSTR( ZAxis,  "ZAxis" )
+   RENUMSTR( CartesianType_end,  "CartesianType_end" )
+RENUMSTR_END()
+
 RENUMSTR_BEGIN( OrbitalType, OrbitalTypeStr )
    RENUMSTR( s,  "s" )
    RENUMSTR( py,  "py" )

--- a/src/base/MathUtilities.h
+++ b/src/base/MathUtilities.h
@@ -68,6 +68,57 @@ template <typename T> T min(T a, T b){
    }
 }
 
+void CalcRotatingMatrix(double matrix[][3], double sita, CartesianType cartesianType){
+   if(cartesianType == XAxis){
+      matrix[0][0] = 1.0;
+      matrix[0][1] = 0.0;
+      matrix[0][2] = 0.0;
+
+      matrix[1][0] = 0.0;
+      matrix[1][1] = cos(sita);
+      matrix[1][2] = sin(sita);
+
+      matrix[2][0] = 0.0;
+      matrix[2][1] = -sin(sita);
+      matrix[2][2] = cos(sita);
+   }
+   else if(cartesianType == YAxis){
+      matrix[0][0] = cos(sita);
+      matrix[0][1] = 0.0;
+      matrix[0][2] = -sin(sita);
+
+      matrix[1][0] = 0.0;
+      matrix[1][1] = 1.0;
+      matrix[1][2] = 0.0;
+
+      matrix[2][0] = sin(sita);
+      matrix[2][1] = 0.0;
+      matrix[2][2] = cos(sita);
+   }
+   else if(cartesianType == ZAxis){
+      matrix[0][0] = cos(sita);
+      matrix[0][1] = sin(sita);
+      matrix[0][2] = 0.0;
+
+      matrix[1][0] = -sin(sita);
+      matrix[1][1] = cos(sita);
+      matrix[1][2] = 0.0;
+
+      matrix[2][0] = 0.0;
+      matrix[2][1] = 0.0;
+      matrix[2][2] = 1.0;
+   }
+   else{
+      stringstream ss;
+      ss << "Error in base::MathUtility::CalcRotatingMatrix: invalid cartesianType \n";
+      throw MolDSException(ss.str());
+   }
+}
+
+
+
+
+
 
 }
 #endif

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -29,7 +29,6 @@ public:
    void OutputConfiguration();
    void CalcPrincipalAxes();
    void Rotate();
-   void Rotate(EularAngle eularAngle);
    void Translate();
 private:
    vector<Atom*>* atomVect;
@@ -39,6 +38,7 @@ private:
    int totalNumberValenceElectrons;
    void CalcInertiaTensor(double** inertiaTensor, double* inertiaTensorOrigin);
    void FreeInertiaTensorMoments(double** inertiaTensor, double* inertiaMoments);
+   void Rotate(EularAngle eularAngle, double* rotatingOrigin, RotatedObjectType rotatedObj);
    void OutputPrincipalAxes(double** inertiaTensor, double* inertiaMoments);
    void OutputInertiaTensorOrigin(double* inertiaTensorOrigin);
    void OutputRotatingConditions(RotatingType rotatingType, double* rotatingOrigin, 
@@ -380,9 +380,9 @@ void Molecule::Rotate(){
    }
    double rotatingOrigin[3] = {this->xyzCOM[0], this->xyzCOM[1], this->xyzCOM[2]};
    if(Parameters::GetInstance()->GetRotatingOrigin() != NULL){
-      rotatingOrigin[0] = Parameters::GetInstance()->GetInertiaTensorOrigin()[0];
-      rotatingOrigin[1] = Parameters::GetInstance()->GetInertiaTensorOrigin()[1];
-      rotatingOrigin[2] = Parameters::GetInstance()->GetInertiaTensorOrigin()[2];
+      rotatingOrigin[0] = Parameters::GetInstance()->GetRotatingOrigin()[0];
+      rotatingOrigin[1] = Parameters::GetInstance()->GetRotatingOrigin()[1];
+      rotatingOrigin[2] = Parameters::GetInstance()->GetRotatingOrigin()[2];
    }
 
    RotatingType rotatingType = Parameters::GetInstance()->GetRotatingType();
@@ -396,16 +396,84 @@ void Molecule::Rotate(){
 
    // rotate
    if(rotatingType == Axis){
+      EularAngle setZAxisEularAngles(rotatingAxis[0], rotatingAxis[1], rotatingAxis[2]);
+      EularAngle angleAroundAxis;
+      angleAroundAxis.SetAlpha(rotatingAngle);
+
+      this->Rotate(setZAxisEularAngles, rotatingOrigin, Frame);
+      this->Rotate(angleAroundAxis, rotatingOrigin, System);
+      this->Rotate(setZAxisEularAngles, rotatingOrigin, System);
    }
    else if(rotatingType == Eular){
-      this->Rotate(rotatingEularAngles);
+      this->Rotate(rotatingEularAngles, rotatingOrigin, System);
    }
-
+   
+   this->OutputConfiguration();
    cout << this->messageDoneRotate;
 }
 
-void Molecule::Rotate(EularAngle eularAngle){
-   // ToDo: rotate
+/***
+ * rotatedObj == System: Molecule is rotated.
+ * rotatedObj == Frame: De Cartesian is rotated.
+ */
+void Molecule::Rotate(EularAngle eularAngle, double* rotatingOrigin, RotatedObjectType rotatedObj){
+
+   double rotatingMatrixAlpha[3][3];
+   double rotatingMatrixBeta[3][3];
+   double rotatingMatrixGamma[3][3];
+   double inv = 1.0;
+   if(rotatedObj == System){
+      inv = -1.0;
+   }
+
+   CalcRotatingMatrix(rotatingMatrixAlpha, inv*eularAngle.GetAlpha(), ZAxis);
+   CalcRotatingMatrix(rotatingMatrixBeta, inv*eularAngle.GetBeta(), YAxis);
+   CalcRotatingMatrix(rotatingMatrixGamma, inv*eularAngle.GetGamma(), ZAxis);
+
+   double temp1[3][3];
+   for(int i=0; i<3; i++){
+      for(int j=0; j<3; j++){
+         temp1[i][j] = 0.0;
+         for(int k=0; k<3; k++){
+            if(rotatedObj == System){
+               temp1[i][j] += rotatingMatrixBeta[i][k] * rotatingMatrixGamma[k][j];
+            }
+            else if(rotatedObj == Frame){
+               temp1[i][j] += rotatingMatrixBeta[i][k] * rotatingMatrixAlpha[k][j];
+            }
+         }
+      }
+   }
+
+   double temp2[3][3];
+   for(int i=0; i<3; i++){
+      for(int j=0; j<3; j++){
+         temp2[i][j] = 0.0;
+         for(int k=0; k<3; k++){
+            if(rotatedObj == System){
+               temp2[i][j] += rotatingMatrixAlpha[i][k] * temp1[k][j];
+            }
+            else if(rotatedObj == Frame){
+               temp2[i][j] += rotatingMatrixGamma[i][k] * temp1[k][j];
+            }
+         }
+      }
+   }
+
+   double rotatedXyz[3];
+   Atom* atom;
+   for(int i=0; i<this->atomVect->size(); i++){
+         atom = (*this->atomVect)[i]; 
+         for(int j=0; j<3; j++){
+            rotatedXyz[j] = 0.0;
+            for(int k=0; k<3; k++){
+               rotatedXyz[j] += temp2[j][k] * (atom->GetXyz()[k] - rotatingOrigin[k]);
+            }
+         }
+         for(int j=0; j<3; j++){
+            atom->GetXyz()[j] = rotatedXyz[j] + rotatingOrigin[j];
+         }
+   }
 }
 
 void Molecule::OutputRotatingConditions(RotatingType rotatingType, double* rotatingOrigin, 

--- a/src/input.in
+++ b/src/input.in
@@ -14,22 +14,27 @@ THEORY
 THEORY_END
 
 INERTIA
-   origin 1.0 2.0 3.0
+   origin 1.0 2.0 30.0
 INERTIA_END
 
 ROTATE
    type axis
    //type eular_angle
-   origin 1.0 2.0 3.0
-   axis 3.0 4.0 5.0
+   //origin 1.0 2.0 3.0
+   origin 0.0 0.0 0.0
+   axis 10.0 0.0 0.0
    angle 30
-   angles 15 25 35
+   //angles 15 25 35
 ROTATE_END
 
 TRANSLATE
    //difference 100.0 200.0 300.0
 TRANSLATE_END
 
+// s
+GEOMETRY
+     S       0.0  0.0  1.0  
+GEOMETRY_END
 
 
 //metane
@@ -61,11 +66,11 @@ TRANSLATE_END
 //GEOMETRY_END
 
 // sh2
-GEOMETRY
-     S       -0.559299    0.471698    0.000000
-     H        0.750701    0.471698    0.000000
-     H       -0.996586    1.706558    0.000000
-GEOMETRY_END
+//GEOMETRY
+//     S       -0.559299    0.471698    0.000000
+//     H        0.750701    0.471698    0.000000
+//     H       -0.996586    1.706558    0.000000
+//GEOMETRY_END
 
 // benzene
 //GEOMETRY