提交: 5986b71a186c91d009337547172cc78ee8da5a30 - MolDS (git)

Commit MetaInfo

• 修订版: 5986b71a186c91d009337547172cc78ee8da5a30 (tree)
• 时间: 2011-01-12 11:01:30
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Logic for core repulsion is added.

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

更改概述

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

差异

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -375,6 +375,7 @@ void InputParser::CalcMolecularBasics(Molecule* molecule){
    molecule->CalcTotalNumberAOs();
    molecule->CalcTotalNumberValenceElectrons();
    molecule->CalcXyzCOM();
+   molecule->CalcTotalCoreRepulsionEnergy();
 
 }
 
@@ -383,6 +384,7 @@ void InputParser::OutputMolecularBasics(Molecule* molecule){
    molecule->OutputTotalNumberAtomsAOsValenceelectrons();
    molecule->OutputConfiguration();
    molecule->OutputXyzCOM();
+   molecule->OutputTotalCoreRepulsionEnergy();
 }
 
 void InputParser::OutputScfConditions(){

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -24,9 +24,11 @@ public:
    void CalcTotalNumberAOs();
    int GetTotalNumberValenceElectrons();
    void CalcTotalNumberValenceElectrons();
+   void CalcTotalCoreRepulsionEnergy();
    void OutputXyzCOM();
    void OutputTotalNumberAtomsAOsValenceelectrons();
    void OutputConfiguration();
+   void OutputTotalCoreRepulsionEnergy();
    void CalcPrincipalAxes();
    void Rotate();
    void Translate();
@@ -38,6 +40,8 @@ private:
    bool wasCalculatedXyzCOM;
    int totalNumberAOs;
    int totalNumberValenceElectrons;
+   double totalCoreRepulsionEnergy;
+   bool wasCalculatedTotalCoreRepulsionEnergy;
    void CalcInertiaTensor(double** inertiaTensor, double* inertiaTensorOrigin);
    void FreeInertiaTensorMoments(double** inertiaTensor, double* inertiaMoments);
    void Rotate(EularAngle eularAngle, double* rotatingOrigin, RotatedObjectType rotatedObj);
@@ -53,6 +57,8 @@ private:
    string messageConfiguration;
    string messageConfigurationTitleAU;
    string messageConfigurationTitleAng;
+   string messageCoreRepulsion;
+   string messageCoreRepulsionTitle;
    string messageCOM;
    string messageCOMTitleAU;
    string messageCOMTitleAng;
@@ -87,12 +93,15 @@ Molecule::Molecule(){
    this->atomVect = new vector<Atom*>;
    this->xyzCOM = MallocerFreer::GetInstance()->MallocDoubleMatrix1d(3);
    this->wasCalculatedXyzCOM = false;
+   this->wasCalculatedTotalCoreRepulsionEnergy = false;
    this->messageTotalNumberAOs = "\tTotal number of valence AOs: ";
    this->messageTotalNumberAtoms = "\tTotal number of atoms: ";
    this->messageTotalNumberValenceElectrons = "\tTotal number of valence electrons: ";
    this->messageConfiguration = "\tMolecular configration:\n";
    this->messageConfigurationTitleAU = "\t\t| i-th | atom type | x [a.u.] | y[a.u.] | z[a.u.] |\n";
    this->messageConfigurationTitleAng = "\t\t| i-th | atom type | x [angst.] | y[angst.] | z[angst.] |\n";
+   this->messageCoreRepulsion = "\tTotal core repulsion energy:\n";
+   this->messageCoreRepulsionTitle = "\t\t| [a.u.] | [eV] |\n";
    this->messageCOM = "\tCenter of Mass:\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";
@@ -194,6 +203,35 @@ int Molecule::GetTotalNumberValenceElectrons(){
    return this->totalNumberValenceElectrons;
 }
 
+void Molecule::CalcTotalCoreRepulsionEnergy(){
+   double energy = 0.0;
+   double distance = 0.0;
+   Atom* atomA = NULL;
+   Atom* atomB = NULL;
+   for(int i=0; i<this->atomVect->size(); i++){
+      atomA = (*this->atomVect)[i];
+      for(int j=i+1; j<this->atomVect->size(); j++){
+         atomB = (*this->atomVect)[j];
+         distance = this->GetDistanceAtoms(i, j);
+         energy += atomA->GetCoreCharge()*atomB->GetCoreCharge()/distance; 
+      }
+   }
+   this->totalCoreRepulsionEnergy = energy;
+   this->wasCalculatedTotalCoreRepulsionEnergy = true;
+}
+
+void Molecule::OutputTotalCoreRepulsionEnergy(){
+   if(!this->wasCalculatedTotalCoreRepulsionEnergy){
+      this->CalcTotalCoreRepulsionEnergy();
+   }
+
+   double eV2AU = Parameters::GetInstance()->GetEV2AU();
+   cout << this->messageCoreRepulsion;
+   cout << this->messageCoreRepulsionTitle;
+   printf("\t\t%e\t%e\n\n",this->totalCoreRepulsionEnergy, this->totalCoreRepulsionEnergy/eV2AU);
+
+}
+   
 void Molecule::CalcTotalNumberValenceElectrons(){
    this->totalNumberValenceElectrons = 0; 
    for(int i=0; i<this->atomVect->size(); i++){

--- a/src/input.in
+++ b/src/input.in
@@ -32,10 +32,10 @@ TRANSLATE
 TRANSLATE_END
 
 // S2
-GEOMETRY
-   S  0.424528    0.741240    0.000000
-   S  -1.353072    0.741240    0.000000
-GEOMETRY_END
+//GEOMETRY
+//   S  0.424528    0.741240    0.000000
+//   S  -1.353072    0.741240    0.000000
+//GEOMETRY_END
 
 // s
 //GEOMETRY
@@ -49,13 +49,13 @@ GEOMETRY_END
 //GEOMETRY_END
 
 //metane
-//GEOMETRY
-//   C  0.647389    0.820131    0.000000
-//   H  1.004043   -0.188679    0.000000
-//   H  1.004062    1.324529    0.873652
-//   H  1.004062    1.324529   -0.873652
-//   H  -0.422611    0.820144    0.000000
-//GEOMETRY_END
+GEOMETRY
+   C  0.647389    0.820131    0.000000
+   H  1.004043   -0.188679    0.000000
+   H  1.004062    1.324529    0.873652
+   H  1.004062    1.324529   -0.873652
+   H  -0.422611    0.820144    0.000000
+GEOMETRY_END
 
 
 // c2