Commit: a20e915b5fa021e81b8d0c4d162bd118846ec811 - MolDS (git)

Commit MetaInfo

• Révision: a20e915b5fa021e81b8d0c4d162bd118846ec811 (tree)
• l'heure: 2012-07-02 19:16:53
• Auteur: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Message de Log

Pm3Pddg::GetDiatomCoreRepulsionSecondDerivative is implemented. #28554

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

Change Summary

• modified: src/pm3/Pm3Pddg.cpp
• modified: src/pm3/Pm3Pddg.h

Modification

--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -181,6 +181,53 @@ double Pm3Pddg::GetDiatomCoreRepulsionFirstDerivative(int atomAIndex,
    return pm3Term + additionalTerm;
 }
 
+// Second derivative of diatomic core repulsion energy.
+// Both derivative are related to the coordinate of atomA.
+double Pm3Pddg::GetDiatomCoreRepulsionSecondDerivative(int atomAIndex,
+                                                       int atomBIndex, 
+                                                       CartesianType axisA1,
+                                                       CartesianType axisA2) const{
+   // PM3 term
+   double pm3Term = Pm3::GetDiatomCoreRepulsionSecondDerivative(atomAIndex, atomBIndex, axisA1, axisA2);
+
+   // pddg additional term, first derivative of eq. (4) in [RCJ_2002]
+   const Atom& atomA = *this->molecule->GetAtom(atomAIndex);
+   const Atom& atomB = *this->molecule->GetAtom(atomBIndex);
+   double distance = this->molecule->GetDistanceAtoms(atomAIndex, atomBIndex);
+   double dCartesian1 = (atomA.GetXyz()[axisA1] - atomB.GetXyz()[axisA1]);
+   double dCartesian2 = (atomA.GetXyz()[axisA2] - atomB.GetXyz()[axisA2]);
+   int na = atomA.GetNumberValenceElectrons();
+   int nb = atomB.GetNumberValenceElectrons();
+   double pddgExponent = -10.0;
+   double tempFirstDeriv = 0.0;
+   double tempSecondDeriv = 0.0;
+   for(int i=0; i<2; i++){
+      double pa = atomA.GetPm3PddgParameterPa(i);
+      double da = atomA.GetPm3PddgParameterDa(i);
+      for(int j=0; j<2; j++){
+         double pb = atomB.GetPm3PddgParameterPa(j);
+         double db = atomB.GetPm3PddgParameterDa(j);
+         tempFirstDeriv +=  this->GetPddgAdditonalDiatomCoreRepulsionTermFirstDerivative(na, pa, da, nb, pb, db, distance);
+         tempSecondDeriv += this->GetPddgAdditonalDiatomCoreRepulsionTermSecondDerivative(na, pa, da, nb, pb, db, distance);
+      }
+   }
+   double preFirstDeriv = 0.0;
+   double preSecondDeriv = 0.0;
+   if(axisA1 != axisA2){
+      preFirstDeriv = -dCartesian1*dCartesian2/pow(distance,3.0);
+      preSecondDeriv = dCartesian1*dCartesian2/pow(distance,2.0);
+   }
+   else{
+      preFirstDeriv = 1.0/distance - dCartesian1*dCartesian1/pow(distance,3.0);
+      preSecondDeriv = pow(dCartesian1/distance,2.0);
+   }
+   preFirstDeriv  /= static_cast<double>(na+nb);
+   preSecondDeriv /= static_cast<double>(na+nb);
+   double additionalTerm = preFirstDeriv*tempFirstDeriv + preSecondDeriv*tempSecondDeriv;
+
+   return pm3Term + additionalTerm;
+}
+
 // see eq. (4) in [RCJ_2002]
 double Pm3Pddg::GetPddgAdditonalDiatomCoreRepulsionTerm(int na, double pa, double da,
                                                         int nb, double pb, double db,

--- a/src/pm3/Pm3Pddg.h
+++ b/src/pm3/Pm3Pddg.h
@@ -34,6 +34,10 @@ protected:
    virtual double GetDiatomCoreRepulsionFirstDerivative(int atomAIndex,
                                                         int atomBIndex, 
                                                         MolDS_base::CartesianType axisA) const;
+   virtual double GetDiatomCoreRepulsionSecondDerivative(int atomAIndex,
+                                                         int atomBIndex, 
+                                                         MolDS_base::CartesianType axisA1,
+                                                         MolDS_base::CartesianType axisA2) const;
 private:
    double GetPddgAdditonalDiatomCoreRepulsionTerm(int na, double pa, double da, 
                                                   int nb, double pb, double db,