Commit: 818049e2c3a2f9916288eb9a3e3fe386515d3b8c - MolDS (git)

Commit MetaInfo

• Revisão: 818049e2c3a2f9916288eb9a3e3fe386515d3b8c (tree)
• Hora: 2014-01-04 17:51:32
• Autor: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Mensagem de Log

Some methods which are not used actually is modified. #32750

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

Mudança Sumário

• modified: src/am1/Am1.cpp
• modified: src/am1/Am1.h
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/mndo/Mndo.cpp
• modified: src/mndo/Mndo.h
• modified: src/pm3/Pm3Pddg.cpp
• modified: src/pm3/Pm3Pddg.h
• modified: src/zindo/ZindoS.cpp

Diff

--- a/src/am1/Am1.cpp
+++ b/src/am1/Am1.cpp
@@ -164,19 +164,16 @@ double Am1::GetDiatomCoreRepulsionEnergy(const Atom& atomA, const Atom& atomB) c
 
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
-double Am1::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                  int indexAtomB, 
-                                                  CartesianType axisA) const{
+double Am1::GetDiatomCoreRepulsion1stDerivative(const Atom& atomA, const Atom& atomB, 
+                                                CartesianType axisA) const{
    // MNDO term
-   double mndoTerms = Mndo::GetDiatomCoreRepulsion1stDerivative(indexAtomA, indexAtomB, axisA);
+   double mndoTerms = Mndo::GetDiatomCoreRepulsion1stDerivative(atomA, atomB, axisA);
 
    // additional term, first derivative of eq. (4) in [S_1989]
    double ang2AU     = Parameters::GetInstance()->GetAngstrom2AU();
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double alphaA     = atomA.GetNddoAlpha(this->theory);
    double alphaB     = atomB.GetNddoAlpha(this->theory);
-   double distance   = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance   = this->molecule->GetDistanceAtoms(atomA, atomB);
    double dCartesian = (atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA]);
    double kA, lA, mA;
    double kB, lB, mB;
@@ -204,20 +201,18 @@ double Am1::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
 
 // Second derivative of diatomic core repulsion energy.
 // Both derivatives are related to the coordinate of atomA.
-double Am1::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                   int indexAtomB, 
-                                                   CartesianType axisA1,
-                                                   CartesianType axisA2) const{
+double Am1::GetDiatomCoreRepulsion2ndDerivative(const Atom& atomA,
+                                                const Atom& atomB, 
+                                                CartesianType axisA1,
+                                                CartesianType axisA2) const{
    // MNDO term
-   double mndoTerm = Mndo::GetDiatomCoreRepulsion2ndDerivative(indexAtomA, indexAtomB, axisA1, axisA2);
+   double mndoTerm = Mndo::GetDiatomCoreRepulsion2ndDerivative(atomA, atomB, axisA1, axisA2);
 
    // additional term, first derivative of eq. (4) in [S_1989]
    double ang2AU     = Parameters::GetInstance()->GetAngstrom2AU();
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double alphaA     = atomA.GetNddoAlpha(this->theory);
    double alphaB     = atomB.GetNddoAlpha(this->theory);
-   double distance   = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance   = this->molecule->GetDistanceAtoms(atomA, atomB);
    double kA, lA, mA;
    double kB, lB, mB;
    double temp1 = 0.0;

--- a/src/am1/Am1.h
+++ b/src/am1/Am1.h
@@ -34,11 +34,11 @@ protected:
    virtual void CalcSCFProperties();
    virtual double GetDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
                                                const MolDS_base_atoms::Atom& atomB) const;
-   virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA) const;
-   virtual double GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion2ndDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA1,
                                                       MolDS_base::CartesianType axisA2) const;
 private:

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -382,12 +382,10 @@ double Cndo2::GetAtomCoreEpcCoulombEnergy(const Atom& atom, const Atom& epc) con
 
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
-double Cndo2::GetDiatomCoreRepulsion1stDerivative(int indexAtomA, int indexAtomB, 
+double Cndo2::GetDiatomCoreRepulsion1stDerivative(const Atom& atomA, const Atom& atomB, 
                                                   CartesianType axisA) const{
    double value=0.0;
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    value = atomA.GetCoreCharge()*atomB.GetCoreCharge();
    value *= (atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA])/distance;
    value *= -1.0/(distance*distance);
@@ -396,8 +394,8 @@ double Cndo2::GetDiatomCoreRepulsion1stDerivative(int indexAtomA, int indexAtomB
 
 // Second derivative of diatomic core repulsion energy.
 // Both derivatives are related to the coordinate of atomA.
-double Cndo2::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                  int indexAtomB, 
+double Cndo2::GetDiatomCoreRepulsion2ndDerivative(const Atom& atomA,
+                                                  const Atom& atomB, 
                                                   CartesianType axisA1,
                                                   CartesianType axisA2) const{
    stringstream ss;
@@ -409,10 +407,10 @@ double Cndo2::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
 void Cndo2::CalcVdWCorrectionEnergy(){
    double value = 0.0;
    for(int i=0; i<this->molecule->GetRealAtomVect().size(); i++){
-      const int indexAtomI = this->molecule->GetRealAtomVect()[i]->GetIndex();
+      const Atom& atomI = *this->molecule->GetAtomVect()[i];
       for(int j=i+1; j<this->molecule->GetRealAtomVect().size(); j++){
-         const int indexAtomJ = this->molecule->GetRealAtomVect()[j]->GetIndex();
-         value += this->GetDiatomVdWCorrectionEnergy(indexAtomI, indexAtomJ);
+         const Atom& atomJ = *this->molecule->GetAtomVect()[j];
+         value += this->GetDiatomVdWCorrectionEnergy(atomI, atomJ);
       }
    }
    this->vdWCorrectionEnergy = value;
@@ -444,10 +442,8 @@ double Cndo2::GetVdwDampingValue2ndDerivative(double vdWDistance, double distanc
 }
 
 // See (2) in [G_2004] ((11) in [G_2006])
-double Cndo2::GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+double Cndo2::GetDiatomVdWCorrectionEnergy(const Atom& atomA, const Atom& atomB) const{
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    double vdWDistance = atomA.GetVdWRadii() + atomB.GetVdWRadii();
    double tmpSum = atomA.GetVdWCoefficient()+atomB.GetVdWCoefficient();
    if(tmpSum<=0e0){return 0e0;}
@@ -460,11 +456,9 @@ double Cndo2::GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const
 
 // First derivative of the vdW correction related to the coordinate of atom A.
 // See (2) in [G_2004] ((11) in [G_2006]).
-double Cndo2::GetDiatomVdWCorrection1stDerivative(int indexAtomA, int indexAtomB, 
+double Cndo2::GetDiatomVdWCorrection1stDerivative(const Atom& atomA, const Atom& atomB, 
                                                   CartesianType axisA) const{
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    double vdWDistance = atomA.GetVdWRadii() + atomB.GetVdWRadii();
    double tmpSum = atomA.GetVdWCoefficient()+atomB.GetVdWCoefficient();
    if(tmpSum<=0e0){return 0e0;}
@@ -485,13 +479,11 @@ double Cndo2::GetDiatomVdWCorrection1stDerivative(int indexAtomA, int indexAtomB
 // Second derivative of the vdW correction.
 // Both derivative sare related to the coordinate of atom A.
 // See (2) in [G_2004] ((11) in [G_2006]).
-double Cndo2::GetDiatomVdWCorrection2ndDerivative(int indexAtomA, 
-                                                     int indexAtomB, 
-                                                     CartesianType axisA1,
-                                                     CartesianType axisA2) const{
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+double Cndo2::GetDiatomVdWCorrection2ndDerivative(const Atom& atomA, 
+                                                  const Atom& atomB, 
+                                                  CartesianType axisA1,
+                                                  CartesianType axisA2) const{
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    double dCartesian1 = atomA.GetXyz()[axisA1] - atomB.GetXyz()[axisA1];
    double dCartesian2 = atomA.GetXyz()[axisA2] - atomB.GetXyz()[axisA2];
    double vdWDistance = atomA.GetVdWRadii() + atomB.GetVdWRadii();

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -136,19 +136,20 @@ protected:
                                                const MolDS_base_atoms::Atom& epc) const;
    virtual double GetDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
                                                const MolDS_base_atoms::Atom& atomB) const;
-   virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA, 
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA) const;
-   virtual double GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion2ndDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA1,
                                                       MolDS_base::CartesianType axisA2) const;
-   virtual double GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const;
-   virtual double GetDiatomVdWCorrection1stDerivative(int indexAtomA, 
-                                                      int indexAtomB, 
+   virtual double GetDiatomVdWCorrectionEnergy(const MolDS_base_atoms::Atom& atomA, 
+                                               const MolDS_base_atoms::Atom& atomB) const;
+   virtual double GetDiatomVdWCorrection1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA) const;
-   virtual double GetDiatomVdWCorrection2ndDerivative(int indexAtomA, 
-                                                      int indexAtomB, 
+   virtual double GetDiatomVdWCorrection2ndDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA1,
                                                       MolDS_base::CartesianType axisA2) const;
    double GetReducedOverlapAOs                      (int na, int nb, double alpha, double beta) const;

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -331,14 +331,11 @@ double Mndo::GetAtomCoreEpcCoulombEnergy(const Atom& atom, const Atom& epc) cons
 
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
-double Mndo::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                 int indexAtomB, 
+double Mndo::GetDiatomCoreRepulsion1stDerivative(const Atom& atomA, const Atom& atomB, 
                                                  CartesianType axisA) const{
    double value =0.0;
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distanceAB = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
-   double twoElecInt = this->twoElecsTwoAtomCores[indexAtomA][indexAtomB][s][s][s][s];
+   double distanceAB = this->molecule->GetDistanceAtoms(atomA, atomB);
+   double twoElecInt = this->twoElecsTwoAtomCores[atomA.GetIndex()][atomB.GetIndex()][s][s][s][s];
    double twoElecInt1stDeriv = this->GetNddoRepulsionIntegral1stDerivative(
                                      atomA, s, s, atomB, s, s, axisA);
    double tmp = this->GetAuxiliaryDiatomCoreRepulsionEnergy(atomA, atomB, distanceAB);
@@ -350,15 +347,13 @@ double Mndo::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
 
 // Second derivative of diatomic core repulsion energy.
 // Both derivatives are related to the coordinate of atomA.
-double Mndo::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                    int indexAtomB, 
-                                                    CartesianType axisA1,
-                                                    CartesianType axisA2) const{
+double Mndo::GetDiatomCoreRepulsion2ndDerivative(const Atom& atomA,
+                                                 const Atom& atomB, 
+                                                 CartesianType axisA1,
+                                                 CartesianType axisA2) const{
    double value =0.0;
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distanceAB = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
-   double twoElecInt = this->twoElecsTwoAtomCores[indexAtomA][indexAtomB][s][s][s][s];
+   double distanceAB = this->molecule->GetDistanceAtoms(atomA, atomB);
+   double twoElecInt = this->twoElecsTwoAtomCores[atomA.GetIndex()][atomB.GetIndex()][s][s][s][s];
    double twoElecInt1stDeriv1 = this->GetNddoRepulsionIntegral1stDerivative(atomA, s, s, 
                                                                             atomB, s, s, 
                                                                             axisA1);
@@ -1595,14 +1590,14 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
          }
 
          // second derivatives of the nuclear repulsions
-         value += this->GetDiatomCoreRepulsion2ndDerivative(indexAtomA, 
-                                                            indexAtomC, 
+         value += this->GetDiatomCoreRepulsion2ndDerivative(atomA, 
+                                                            atomC, 
                                                             static_cast<CartesianType>(axisA1), 
                                                             static_cast<CartesianType>(axisA2));
          // second derivatives of the van der waals corrections
          if(Parameters::GetInstance()->RequiresVdWSCF()){
-            value += this->GetDiatomVdWCorrection2ndDerivative(indexAtomA, 
-                                                               indexAtomC, 
+            value += this->GetDiatomVdWCorrection2ndDerivative(atomA, 
+                                                               atomC, 
                                                                static_cast<CartesianType>(axisA1), 
                                                                static_cast<CartesianType>(axisA2));
          }
@@ -1759,14 +1754,14 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
    }
 
    // second derivatives of the nuclear repulsions
-   value -= this->GetDiatomCoreRepulsion2ndDerivative(indexAtomA, 
-                                                      indexAtomB, 
+   value -= this->GetDiatomCoreRepulsion2ndDerivative(atomA, 
+                                                      atomB, 
                                                       static_cast<CartesianType>(axisA), 
                                                       static_cast<CartesianType>(axisB));
    // second derivatives of the van der waals corrections
    if(Parameters::GetInstance()->RequiresVdWSCF()){
-      value -= this->GetDiatomVdWCorrection2ndDerivative(indexAtomA, 
-                                                         indexAtomB, 
+      value -= this->GetDiatomVdWCorrection2ndDerivative(atomA, 
+                                                         atomB, 
                                                          static_cast<CartesianType>(axisA), 
                                                          static_cast<CartesianType>(axisB));
    }
@@ -2723,10 +2718,10 @@ void Mndo::CalcForce(const vector<int>& elecStates){
                double coreRepulsion[CartesianType_end] = {0.0,0.0,0.0};
                for(int i=0; i<CartesianType_end; i++){
                   coreRepulsion[i] += this->GetDiatomCoreRepulsion1stDerivative(
-                                            a, b, (CartesianType)i);
+                                            atomA, atomB, (CartesianType)i);
                   if(Parameters::GetInstance()->RequiresVdWSCF()){
                      coreRepulsion[i] += this->GetDiatomVdWCorrection1stDerivative(
-                                               a, b, (CartesianType)i);
+                                               atomA, atomB, (CartesianType)i);
                   }
                }  
                // electron core attraction part (ground state)

--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -56,11 +56,11 @@ protected:
                                                const MolDS_base_atoms::Atom& epc) const;
    virtual double GetDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
                                                const MolDS_base_atoms::Atom& atomB) const;
-   virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA) const;
-   virtual double GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion2ndDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA1,
                                                       MolDS_base::CartesianType axisA2) const;
    virtual double GetFockDiagElement(const MolDS_base_atoms::Atom& atomA, 

--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -163,18 +163,15 @@ double Pm3Pddg::GetDiatomCoreRepulsionEnergy(const Atom& atomA, const Atom& atom
 
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
-double Pm3Pddg::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                    int indexAtomB, 
+double Pm3Pddg::GetDiatomCoreRepulsion1stDerivative(const Atom& atomA, const Atom& atomB, 
                                                     CartesianType axisA) const{
    // PM3 term
-   double pm3Term = Pm3::GetDiatomCoreRepulsion1stDerivative(indexAtomA, indexAtomB, axisA);
+   double pm3Term = Pm3::GetDiatomCoreRepulsion1stDerivative(atomA, atomB, axisA);
 
    // pddg additional term, first derivative of eq. (4) in [RCJ_2002]
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int na = atomA.GetNumberValenceElectrons();
    int nb = atomB.GetNumberValenceElectrons();
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    double temp = 0.0;
    for(int i=0; i<2; i++){
       double pa = atomA.GetPm3PddgParameterPa(i);
@@ -193,17 +190,15 @@ double Pm3Pddg::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
 
 // Second derivative of diatomic core repulsion energy.
 // Both derivative are related to the coordinate of atomA.
-double Pm3Pddg::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                       int indexAtomB, 
-                                                       CartesianType axisA1,
-                                                       CartesianType axisA2) const{
+double Pm3Pddg::GetDiatomCoreRepulsion2ndDerivative(const Atom& atomA,
+                                                    const Atom& atomB, 
+                                                    CartesianType axisA1,
+                                                    CartesianType axisA2) const{
    // PM3 term
-   double pm3Term = Pm3::GetDiatomCoreRepulsion2ndDerivative(indexAtomA, indexAtomB, axisA1, axisA2);
+   double pm3Term = Pm3::GetDiatomCoreRepulsion2ndDerivative(atomA, atomB, axisA1, axisA2);
 
    // pddg additional term, first derivative of eq. (4) in [RCJ_2002]
-   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
-   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
-   double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
+   double distance = this->molecule->GetDistanceAtoms(atomA, atomB);
    double dCartesian1 = (atomA.GetXyz()[axisA1] - atomB.GetXyz()[axisA1]);
    double dCartesian2 = (atomA.GetXyz()[axisA2] - atomB.GetXyz()[axisA2]);
    int na = atomA.GetNumberValenceElectrons();

--- a/src/pm3/Pm3Pddg.h
+++ b/src/pm3/Pm3Pddg.h
@@ -32,11 +32,11 @@ protected:
    virtual void SetEnableAtomTypes();
    virtual double GetDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
                                                const MolDS_base_atoms::Atom& atomB) const;
-   virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA) const;
-   virtual double GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
-                                                      int indexAtomB, 
+   virtual double GetDiatomCoreRepulsion2ndDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                      const MolDS_base_atoms::Atom& atomB, 
                                                       MolDS_base::CartesianType axisA1,
                                                       MolDS_base::CartesianType axisA2) const;
 private:

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -3717,9 +3717,9 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                double forceElecCoreAttPart[CartesianType_end] = {0.0,0.0,0.0};
                for(int i=0; i<CartesianType_end; i++){
                   // core repulsion part (ground state)
-                  coreRepulsion[i] = this->GetDiatomCoreRepulsion1stDerivative(a, b, static_cast<CartesianType>(i));
+                  coreRepulsion[i] = this->GetDiatomCoreRepulsion1stDerivative(atomA, atomB, static_cast<CartesianType>(i));
                   if(Parameters::GetInstance()->RequiresVdWSCF()){
-                     coreRepulsion[i] += this->GetDiatomVdWCorrection1stDerivative(a, b, static_cast<CartesianType>(i));
+                     coreRepulsion[i] += this->GetDiatomVdWCorrection1stDerivative(atomA, atomB, static_cast<CartesianType>(i));
                   }
                   // electron core attraction part (ground state)                     
                   forceElecCoreAttPart[i] = ( atomA.GetCoreCharge()*atomicElectronPopulation[b]
@@ -3867,10 +3867,10 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
 
                   // Calculation of core repusion force
                   coreRepulsion += this->GetDiatomCoreRepulsion1stDerivative(
-                                         a, b, (CartesianType)i);
+                                         atomA, atomB, (CartesianType)i);
                   if(Parameters::GetInstance()->RequiresVdWSCF()){
                      coreRepulsion += this->GetDiatomVdWCorrection1stDerivative(
-                                            a, b, (CartesianType)i);
+                                            atomA, atomB, (CartesianType)i);
                   }
 
                   // Calculate force arise from electronic part.