Commit: 592284608ac794c52cea232f74f41bbfa41ae1cf - MolDS (git)

Commit MetaInfo

• Revisión: 592284608ac794c52cea232f74f41bbfa41ae1cf (tree)
• Tiempo: 2013-05-22 07:32:17
• Autor: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

ZindoS::CalcForce for excited states is paritally implemented although not completed yet. #31221

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

Cambiar Resumen

• modified: src/zindo/ZindoS.cpp
• modified: src/zindo/ZindoS.h

Diferencia incremental

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -3937,6 +3937,31 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                // excited state force
                for(int n=0; n<elecStates.size(); n++){
                   if(elecStates[n]<=0){continue;}
+                  // static part
+                  double forceExcitedStaticPart[CartesianType_end] = {0.0,0.0,0.0};
+                  this->CalcForceExcitedStaticPart(forceExcitedStaticPart,
+                                                   n,
+                                                   a,
+                                                   b,
+                                                   diatomicTwoElecTwoCore1stDerivs);
+                  // sum up contributions from static part (excited state)
+#pragma omp critical
+                  {
+                     for(int i=0; i<CartesianType_end; i++){
+                        this->matrixForce[n][b][i] += forceExcitedStaticPart[i];
+                        this->matrixForce[n][a][i] -= forceExcitedStaticPart[i];
+                     }
+                  }
+                  
+                  // response part
+                  // electron core attraction part (excited state)
+                  double forceExcitedElecCoreAttPart[CartesianType_end]={0.0,0.0,0.0};
+                  this->CalcForceExcitedElecCoreAttractionPart(
+                                             forceExcitedElecCoreAttPart,
+                                             n,
+                                             a,
+                                             b,
+                                             diatomicTwoElecTwoCore1stDerivs);
                }
             } // end of for(int b)
          }    // end of for(int a)
@@ -3960,7 +3985,8 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
    }
   
    /*
-   // Calculate force. First derivative of overlapAOs integral is
+   // Calculate force (on the ground state only). 
+   // First derivative of overlapAOs integral is
    // calculated with GTO expansion technique.
    stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
@@ -4040,5 +4066,52 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
    }
    */
 }
+
+void ZindoS::CalcForceExcitedStaticPart(double* force, 
+                                      int elecStateIndex,
+                                      int indexAtomA, 
+                                      int indexAtomB,
+                                      double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
+   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   int firstAOIndexA = atomA.GetFirstAOIndex();
+   int firstAOIndexB = atomB.GetFirstAOIndex();
+   int lastAOIndexA  = atomA.GetLastAOIndex();
+   int lastAOIndexB  = atomB.GetLastAOIndex();
+   for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
+      for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
+         for(int i=0; i<CartesianType_end; i++){
+            double temp= 2.0*this->etaMatrixForce[elecStateIndex][mu][mu]
+                            *this->etaMatrixForce[elecStateIndex][lambda][lambda]
+                        -1.0*this->etaMatrixForce[elecStateIndex][mu][lambda]
+                            *this->etaMatrixForce[elecStateIndex][mu][lambda];
+            force[i] += temp
+                       *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
+                                                       [lambda-firstAOIndexB]
+                                                       [i];
+         }
+      }
+   }
+}
+
+void ZindoS::CalcForceExcitedElecCoreAttractionPart(double* force, 
+                                                  int elecStateIndex,
+                                                  int indexAtomA, 
+                                                  int indexAtomB,
+                                                  double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
+   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   int firstAOIndexA = atomA.GetFirstAOIndex();
+   int lastAOIndexA  = atomA.GetLastAOIndex();
+   for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
+      for(int i=0; i<CartesianType_end; i++){
+         force[i] += this->zMatrixForce[elecStateIndex][mu][mu]
+                    *atomB.GetCoreCharge()
+                    *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA][s][i];
+      }
+   }
+}
+
+
 }
 

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -237,6 +237,16 @@ private:
    void CalcDiatomicTwoElecTwoCore1stDerivatives(double*** matrix, 
                                                  int indexAtomA, 
                                                  int indexAtomB) const;
+   void CalcForceExcitedStaticPart(double* force, 
+                                   int elecStateIndex,
+                                   int indexAtomA,
+                                   int indexAtomB,
+                                   double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
+   void CalcForceExcitedElecCoreAttractionPart(double* force, 
+                                               int elecStateIndex,
+                                               int indexAtomA,
+                                               int indexAtomB,
+                                               double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
    void CalcFreeExcitonEnergies(double** freeExcitonEnergiesCIS, 
                                 const MolDS_base::Molecule& molecule, 
                                 double const* energiesMO,