提交: 95987bb7e7082972fe73624d76c1815b2c648d78 - MolDS (git)

Commit MetaInfo

• 修订版: 95987bb7e7082972fe73624d76c1815b2c648d78 (tree)
• 时间: 2011-11-25 22:56:53
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Atom::GetAtomicBasisValue is implemented. #26391

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

更改概述

• modified: src/base/EularAngle.cpp
• modified: src/base/atoms/Atom.cpp
• modified: src/base/atoms/Atom.h

差异

--- a/src/base/EularAngle.cpp
+++ b/src/base/EularAngle.cpp
@@ -11,8 +11,9 @@ namespace MolDS_base{
 
 EularAngle::EularAngle(){
    this->SetMessage();
-   this->alpha = 0.0;
-   this->beta = 0.0;
+   // e the [BFB_1997] for defenitions of alpha, beta, gamma;
+   this->alpha = 0.0; // (= "phi" in P25 in J. A. Pople book)
+   this->beta = 0.0;  // (= "theta" in P25 in J. A. Pople book)
    this->gamma = 0.0;
 }
 

--- a/src/base/atoms/Atom.cpp
+++ b/src/base/atoms/Atom.cpp
@@ -8,7 +8,9 @@
 #include<stdexcept>
 #include"../MolDSException.h"
 #include"../Enums.h"
+#include"../MathUtilities.h"
 #include"../MallocerFreer.h"
+#include"../EularAngle.h"
 #include"Atom.h"
 using namespace std;
 using namespace MolDS_base;
@@ -51,6 +53,12 @@ void Atom::SetMessages(){
    this->errorMessageOrbitalType = "\torbital type = ";
    this->errorMessageShellType = "\tshell type = ";
    this->errorMessageTheoryType = "\tTheory = ";
+   this->errorMessageNumberValences = "\tnumber of valences = ";
+   this->errorMessageValenceIndex = "\tvalenceIndex = ";
+   this->errorMessageGetAtomicBasisValueBadValenceIndex 
+      = "Error in molds_atoms::Atom::GetAtomicBasisValue: Bad valenceIndex is set.\n";
+   this->errorMessageGetRealAnuglarPartAOBadValence 
+      = "Error in molds_atoms::Atom::GetRealAnuglarPartAO: Bad valence orbital is set.\n";
    this->errorMessageEffectivPrincipalQuantumNumber = 
       "Error in base::Atom::GetEffectivePrincipalQuantumNumber: invalid shelltype.\n";
    this->errorMessageZindoCoreIntegral = "Error in base_atoms::Atom::GetZindoCoreINtegral: Invalid orbitalType.\n";
@@ -129,6 +137,85 @@ vector<OrbitalType> Atom::GetValence(){
    return this->valence;
 }
 
+double Atom::GetAtomicBasisValue(double x, 
+                                 double y, 
+                                 double z, 
+                                 int valenceIndex, 
+                                 TheoryType theory){
+   if(this->valence.size()<=valenceIndex){
+      stringstream ss;
+      ss << this->errorMessageGetAtomicBasisValueBadValenceIndex;
+      ss << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+      ss << this->errorMessageNumberValences << this->valence.size() << endl;
+      ss << this->errorMessageValenceIndex << valenceIndex << endl;
+      throw MolDSException(ss.str());
+   }
+   double dx = x - this->xyz[XAxis];
+   double dy = y - this->xyz[YAxis];
+   double dz = z - this->xyz[ZAxis];
+   double dr = sqrt( pow(dx,2.0) + pow(dy,2.0) + pow(dz,2.0) );
+   EularAngle eularAngle(dx, dy, dz);
+   double angularPart = this->GetRealAnuglarPartAO(eularAngle.GetBeta(),
+                                                   eularAngle.GetAlpha(),
+                                                   this->valence[valenceIndex]);
+   double orbitalExponent = this->GetOrbitalExponent(this->valenceShellType,
+                                                     this->valence[valenceIndex],
+                                                     theory);
+   double radialPart = this->GetRadialPartAO(dr, orbitalExponent, this->valenceShellType);
+   return angularPart*radialPart;
+}
+
+// See (1.74) & (1.72) in J. A. Pople book.
+double Atom::GetRadialPartAO(double dr, double orbitalExponent, MolDS_base::ShellType shell){
+   int principalQuantumNumber = (int)shell + 1;
+   double temp1 = pow(2.0*orbitalExponent,(double)principalQuantumNumber+0.5);
+   double temp2 = pow(Factorial(2*principalQuantumNumber),-0.5);
+   return temp1*temp2*pow(dr,principalQuantumNumber-1)*exp(-1.0*orbitalExponent*dr);
+}
+
+// See Table 1 in [BFB_1997] or Table 1.2 in J. A. Pople book.
+// See Table 1 in [BFB_1997] or p25 in J. A. Pople book for defenitions of theta and phi.
+double Atom::GetRealAnuglarPartAO(double theta, double phi, OrbitalType orbital){
+   double value=0.0;
+   switch(orbital){
+      case s:
+         value = pow(4.0*M_PI,-0.5);
+         break;
+      case py:
+         value = pow(3.0/(4.0*M_PI),0.5)*sin(theta)*sin(phi);
+         break;
+      case pz:
+         value = pow(3.0/(4.0*M_PI),0.5)*cos(theta);
+         break;
+      case px:
+         value = pow(3.0/(4.0*M_PI),0.5)*sin(theta)*cos(phi);
+         break;
+      case dxy:
+         value = pow(15.0/(16.0*M_PI),0.5)*pow(sin(theta),2.0)*sin(2.0*phi);
+         break;
+      case dyz:
+         value = pow(15.0/(16.0*M_PI),0.5)*sin(2.0*theta)*sin(phi);
+         break;
+      case dzz:
+         value = pow(5.0/(16.0*M_PI),0.5)*(3.0*pow(cos(theta),2.0) - 1.0);
+         break;
+      case dzx:
+         value = pow(15.0/(16.0*M_PI),0.5)*sin(2.0*theta)*cos(phi);
+         break;
+      case dxxyy:
+         value = pow(15.0/(16.0*M_PI),0.5)*pow(sin(theta),2.0)*cos(2.0*phi);
+         break;
+      default:
+         stringstream ss;
+         ss << this->errorMessageGetRealAnuglarPartAOBadValence;
+         ss << this->errorMessageAtomType << AtomTypeStr(this->atomType) << endl;
+         ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbital) << endl;
+         throw MolDSException(ss.str());
+   }
+   return value;
+}
+
+
 double Atom::GetBondingParameter(TheoryType theory, OrbitalType orbital){
 
    double value = 0.0;

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -14,6 +14,11 @@ public:
    double* GetPxyz();
    void SetPxyz(double px, double py, double pz);
    std::vector<MolDS_base::OrbitalType> GetValence();
+   double GetAtomicBasisValue(double x, 
+                              double y, 
+                              double z, 
+                              int valenceIndex,
+                              MolDS_base::TheoryType theory);
    double GetBondingParameter();
    double GetBondingParameter(MolDS_base::TheoryType theory, 
                               MolDS_base::OrbitalType orbital);
@@ -167,9 +172,10 @@ protected:
    double pm3ParameterL[4];// Table II in ref. [S_1989].
    double pm3ParameterM[4];// Table II in ref. [S_1989].
 private:
-   void SetMessages();
    std::string errorMessageIonPot;
    std::string errorMessageAtomType;
+   std::string errorMessageNumberValences;
+   std::string errorMessageValenceIndex;
    std::string errorMessageOrbitalType;
    std::string errorMessageOrbitalExponent;
    std::string errorMessageShellType;
@@ -180,6 +186,8 @@ private:
    std::string errorMessageMndoCoreIntegral;
    std::string errorMessageAm1CoreIntegral;
    std::string errorMessagePm3CoreIntegral;
+   std::string errorMessageGetAtomicBasisValueBadValenceIndex;
+   std::string errorMessageGetRealAnuglarPartAOBadValence;
    std::string errorMessageGetOrbitalExponentBadTheory;
    std::string errorMessageTheoryType;
    std::string errorMessageGetBondingParameterBadTheoryBadOrbital;
@@ -205,6 +213,9 @@ private:
    std::string errorMessageGetNddoGpp2BadTheory;
    std::string errorMessageGetNddoHspBadTheory;
    std::string errorMessageGetNddoHppBadTheory;
+   void SetMessages();
+   double GetRealAnuglarPartAO(double theta, double phi, MolDS_base::OrbitalType orbital);
+   double GetRadialPartAO(double dr, double orbitalExponent, MolDS_base::ShellType shell);
    int GetEffectivePrincipalQuantumNumber(MolDS_base::ShellType shellType); // Table 1.4 in J. A. Pople book
    double GetZindoJss();  // Part of Eq. (13) in [BZ_1979]
    double GetZindoJsp();  // Part of Eq. (13) in [BZ_1979]