Commit: b1162b09154403db21876430989860ecb4b9d702 - MolDS (git)

Commit MetaInfo

• Revision: b1162b09154403db21876430989860ecb4b9d702 (tree)
• Zeit: 2011-09-30 22:29:35
• Autor: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Center of core is added.

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

Ändern Zusammenfassung

• modified: src/base/InputParser.h
• modified: src/base/Molecule.h
• modified: src/input.in

Diff

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -589,6 +589,7 @@ void InputParser::CalcMolecularBasics(Molecule* molecule){
    molecule->CalcTotalNumberAOs();
    molecule->CalcTotalNumberValenceElectrons();
    molecule->CalcXyzCOM();
+   molecule->CalcXyzCOC();
    molecule->CalcTotalCoreRepulsionEnergy();
 
 }
@@ -638,6 +639,7 @@ void InputParser::OutputMolecularBasics(Molecule* molecule){
    molecule->OutputTotalNumberAtomsAOsValenceelectrons();
    molecule->OutputConfiguration();
    molecule->OutputXyzCOM();
+   molecule->OutputXyzCOC();
    molecule->OutputTotalCoreRepulsionEnergy();
 }
 

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -19,7 +19,9 @@ public:
    ~Molecule();
    vector<Atom*>* GetAtomVect(); 
    double* GetXyzCOM();
+   double* GetXyzCOC();
    void CalcXyzCOM();
+   void CalcXyzCOC();
    int GetTotalNumberAOs();
    void CalcTotalNumberAOs();
    int GetTotalNumberValenceElectrons();
@@ -28,6 +30,7 @@ public:
    void CalcTotalCoreRepulsionEnergy();
    double GetTotalCoreRepulsionEnergy();
    void OutputXyzCOM();
+   void OutputXyzCOC();
    void OutputTotalNumberAtomsAOsValenceelectrons();
    void OutputConfiguration();
    void OutputTotalCoreRepulsionEnergy();
@@ -38,8 +41,10 @@ public:
    double GetDistanceAtoms(Atom* atomA, Atom* atomB);
 private:
    vector<Atom*>* atomVect;
-   double* xyzCOM;
+   double* xyzCOM; // x, y, z coordinates of Center of Mass;
+   double* xyzCOC; // x, y, z coordinates of Center of Core;
    bool wasCalculatedXyzCOM;
+   bool wasCalculatedXyzCOC;
    int totalNumberAOs;
    int totalNumberValenceElectrons;
    double totalCoreRepulsionEnergy;
@@ -62,6 +67,7 @@ private:
    string messageCoreRepulsion;
    string messageCoreRepulsionTitle;
    string messageCOM;
+   string messageCOC;
    string messageCOMTitleAU;
    string messageCOMTitleAng;
    string messageStartPrincipalAxes;
@@ -94,7 +100,9 @@ private:
 Molecule::Molecule(){
    this->atomVect = new vector<Atom*>;
    this->xyzCOM = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
+   this->xyzCOC = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
    this->wasCalculatedXyzCOM = false;
+   this->wasCalculatedXyzCOC = false;
    this->wasCalculatedTotalCoreRepulsionEnergy = false;
    this->messageTotalNumberAOs = "\tTotal number of valence AOs: ";
    this->messageTotalNumberAtoms = "\tTotal number of atoms: ";
@@ -105,6 +113,7 @@ Molecule::Molecule(){
    this->messageCoreRepulsion = "\tTotal core repulsion energy:\n";
    this->messageCoreRepulsionTitle = "\t\t| [a.u.] | [eV] |\n";
    this->messageCOM = "\tCenter of Mass:\n";
+   this->messageCOC = "\tCenter of Core:\n";
    this->messageCOMTitleAU = "\t\t| x [a.u.] | y[a.u.] | z[a.u.] |\n";
    this->messageCOMTitleAng = "\t\t| x [angst.] | y[angst.] | z[angst.] |\n";
    this->messageStartPrincipalAxes = "**********  START: Principal Axes of Inertia  **********\n";
@@ -148,6 +157,10 @@ Molecule::~Molecule(){
       MallocerFreer::GetInstance()->FreeDoubleMatrix1d(&this->xyzCOM);
       //cout << "xyzCOM deleted\n";
    }
+   if(this->xyzCOC != NULL){
+      MallocerFreer::GetInstance()->FreeDoubleMatrix1d(&this->xyzCOC);
+      //cout << "xyzCOC deleted\n";
+   }
 }
 
 vector<Atom*>* Molecule::GetAtomVect(){
@@ -161,11 +174,18 @@ double* Molecule::GetXyzCOM(){
    return this->xyzCOM;
 }
 
+double* Molecule::GetXyzCOC(){
+   if(!this->wasCalculatedXyzCOC){
+      this->CalcXyzCOC();
+   }
+   return this->xyzCOC;
+}
+
 void Molecule::CalcXyzCOM(){
-   double totalAtomicMass;
+   double totalAtomicMass = 0.0;
    Atom* atom;
    double* atomicXyz;
-   double atomicMass;
+   double atomicMass = 0.0;
 
    for(int j=0; j<3; j++){
       this->xyzCOM[j] = 0.0;
@@ -186,6 +206,31 @@ void Molecule::CalcXyzCOM(){
    this->wasCalculatedXyzCOM = true;
 }
 
+void Molecule::CalcXyzCOC(){
+   double totalCoreMass = 0.0;
+   Atom* atom;
+   double* atomicXyz;
+   double coreMass = 0.0;
+
+   for(int j=0; j<3; j++){
+      this->xyzCOC[j] = 0.0;
+   }
+      
+   for(int i=0; i<this->atomVect->size(); i++){
+      atom = (*this->atomVect)[i]; 
+      atomicXyz = atom->GetXyz();
+      coreMass = atom->GetAtomicMass() - (double)atom->GetNumberValenceElectrons();
+      totalCoreMass += coreMass;
+      for(int j=0; j<3; j++){
+         this->xyzCOC[j] += atomicXyz[j] * coreMass;
+      }
+   }
+   for(int i=0; i<3; i++){
+      this->xyzCOC[i]/=totalCoreMass;
+   }
+   this->wasCalculatedXyzCOC = true;
+}
+
 int Molecule::GetTotalNumberAOs(){
    return this->totalNumberAOs;
 }
@@ -296,6 +341,23 @@ void Molecule::OutputXyzCOM(){
 
 }
 
+void Molecule::OutputXyzCOC(){
+   double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();
+   cout << this->messageCOC;
+   cout << this->messageCOMTitleAng;
+   printf("\t\t%e\t%e\t%e\n",this->xyzCOC[0]/ang2AU,
+                             this->xyzCOC[1]/ang2AU,
+                             this->xyzCOC[2]/ang2AU);
+   cout << "\n";
+
+   cout << this->messageCOMTitleAU;
+   printf("\t\t%e\t%e\t%e\n",this->xyzCOC[0],
+                             this->xyzCOC[1],
+                             this->xyzCOC[2]);
+   cout << "\n";
+
+}
+
 void Molecule::OutputTotalNumberAtomsAOsValenceelectrons(){
    cout << this->messageTotalNumberAtoms << this->atomVect->size() << "\n";
    cout << this->messageTotalNumberAOs << this->totalNumberAOs << "\n";
@@ -603,9 +665,12 @@ void Molecule::Translate(){
    
    this->wasCalculatedXyzCOM = false;
    this->CalcXyzCOM();
+   this->wasCalculatedXyzCOC = false;
+   this->CalcXyzCOC();
 
    this->OutputConfiguration();
    this->OutputXyzCOM();
+   this->OutputXyzCOC();
 
    cout << this->messageDoneTranslate;
 }

--- a/src/input.in
+++ b/src/input.in
@@ -189,6 +189,13 @@ GEOMETRY_END
 //   Li  1.994960      0.485175      0.000000
 //GEOMETRY_END
 
+// LiH 3
+//GEOMETRY
+//   H   -1.14750   0.0      0.0
+//   Li  0.382500   0.0      0.0
+//GEOMETRY_END
+
+
 // Li2
 //GEOMETRY
 //   Li  0.0   0.0   1.230000