OSDN -- Desenvolver e Download de Software Open Source
Estado tradução de Português translation status for pt
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
No Tags

Frequently used words (click to add to your profile)

javac++androidlinuxc#windowsobjective-ccocoa誰得qtpythonphprubygameguibathyscaphec計画中(planning stage)翻訳omegatframeworktwitterdomtestvb.netdirectxゲームエンジンbtronarduinopreviewer

Commit MetaInfo

Revisão209b733d771777da7f71f3df2377b4e2d028a8f6 (tree)
Hora2013-05-22 10:53:12
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Mensagem de Log

ZindoS::Calcforce is modified for excited states, although not completed yet. #31221

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

Mudança Sumário

Diff

--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -1369,7 +1369,7 @@ void InputParser::ValidateMdConditions(const Molecule& molecule) const{
13691369 int targetStateIndex = Parameters::GetInstance()->GetElectronicStateIndexMD();
13701370 TheoryType theory = Parameters::GetInstance()->GetCurrentTheory();
13711371 // Validate theory
1372- if(theory == CNDO2 || theory == INDO || (theory == ZINDOS && groundStateIndex < targetStateIndex)){
1372+ if(theory == CNDO2 || theory == INDO ){
13731373 stringstream ss;
13741374 ss << this->errorMessageNonValidExcitedStatesMD;
13751375 ss << this->errorMessageElecState << targetStateIndex << endl;
--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -2196,36 +2196,6 @@ void Mndo::CalcForceExcitedElecCoreAttractionPart(double* force,
21962196 }
21972197 }
21982198
2199-void Mndo::CalcForceExcitedOverlapAOsPart(double* force,
2200- int elecStateIndex,
2201- int indexAtomA,
2202- int indexAtomB,
2203- double const* const* const* diatomicOverlapAOs1stDerivs) const{
2204- const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
2205- const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
2206- int firstAOIndexA = atomA.GetFirstAOIndex();
2207- int firstAOIndexB = atomB.GetFirstAOIndex();
2208- int lastAOIndexA = atomA.GetLastAOIndex();
2209- int lastAOIndexB = atomB.GetLastAOIndex();
2210- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2211- for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
2212- double bondParameter = atomA.GetBondingParameter(
2213- this->theory,
2214- atomA.GetValence(mu-firstAOIndexA))
2215- +atomB.GetBondingParameter(
2216- this->theory,
2217- atomB.GetValence(nu-firstAOIndexB));
2218- bondParameter *= 0.5;
2219- for(int i=0; i<CartesianType_end; i++){
2220- force[i] += -1.0
2221- *this->zMatrixForce[elecStateIndex][mu][nu]
2222- *bondParameter
2223- *diatomicOverlapAOs1stDerivs[mu-firstAOIndexA][nu-firstAOIndexB][i];
2224- }
2225- }
2226- }
2227-}
2228-
22292199 void Mndo::CalcForceExcitedTwoElecPart(double* force,
22302200 int elecStateIndex,
22312201 int indexAtomA,
@@ -2358,7 +2328,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
23582328 }
23592329
23602330 // response part
2361- // electron core attraction part (excited state)
2331+ // electron core attraction part (excited states)
23622332 double forceExcitedElecCoreAttPart[CartesianType_end]={0.0,0.0,0.0};
23632333 this->CalcForceExcitedElecCoreAttractionPart(
23642334 forceExcitedElecCoreAttPart,
@@ -2366,14 +2336,14 @@ void Mndo::CalcForce(const vector<int>& elecStates){
23662336 a,
23672337 b,
23682338 diatomicTwoElecTwoCore1stDerivs);
2369- // overlapAOs part (excited state)
2339+ // overlapAOs part (excited states)
23702340 double forceExcitedOverlapAOsPart[CartesianType_end] = {0.0,0.0,0.0};
23712341 this->CalcForceExcitedOverlapAOsPart(forceExcitedOverlapAOsPart,
23722342 n,
23732343 a,
23742344 b,
23752345 diatomicOverlapAOs1stDerivs);
2376- // two electron part (ground state)
2346+ // two electron part (excited states)
23772347 double forceExcitedTwoElecPart[CartesianType_end] = {0.0,0.0,0.0};
23782348 this->CalcForceExcitedTwoElecPart(forceExcitedTwoElecPart,
23792349 n,
@@ -2434,6 +2404,190 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs, dou
24342404 CartesianType_end);
24352405 }
24362406
2407+// see (18) in [PT_1997]
2408+double Mndo::GetSmallQElement(int moI,
2409+ int moP,
2410+ double const* const* xiOcc,
2411+ double const* const* xiVir,
2412+ double const* const* eta) const{
2413+ double value = 0.0;
2414+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
2415+ bool isMoPOcc = moP<numberOcc ? true : false;
2416+
2417+ for(int A=0; A<molecule->GetNumberAtoms(); A++){
2418+ const Atom& atomA = *molecule->GetAtom(A);
2419+ int firstAOIndexA = atomA.GetFirstAOIndex();
2420+ int lastAOIndexA = atomA.GetLastAOIndex();
2421+
2422+ for(int B=A; B<molecule->GetNumberAtoms(); B++){
2423+ const Atom& atomB = *molecule->GetAtom(B);
2424+ int firstAOIndexB = atomB.GetFirstAOIndex();
2425+ int lastAOIndexB = atomB.GetLastAOIndex();
2426+
2427+ if(A!=B){
2428+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2429+ for(int nu=mu; nu<=lastAOIndexA; nu++){
2430+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2431+ for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
2432+ double twoElecInt = 0.0;
2433+ twoElecInt = this->twoElecTwoCore[A]
2434+ [B]
2435+ [mu-firstAOIndexA]
2436+ [nu-firstAOIndexA]
2437+ [lambda-firstAOIndexB]
2438+ [sigma-firstAOIndexB];
2439+ double temp = 0.0;
2440+ if(isMoPOcc){
2441+ int p = numberOcc - (moP+1);
2442+ temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
2443+ -1.0*xiOcc[p][lambda]*eta[nu][sigma]
2444+ -1.0*xiOcc[p][sigma]*eta[nu][lambda];
2445+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
2446+ temp = 4.0*xiOcc[p][sigma]*eta[mu][nu]
2447+ -1.0*xiOcc[p][mu]*eta[sigma][nu]
2448+ -1.0*xiOcc[p][nu]*eta[sigma][mu];
2449+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
2450+ }
2451+ else{
2452+ int p = moP - numberOcc;
2453+ temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
2454+ -1.0*xiVir[p][lambda]*eta[sigma][nu]
2455+ -1.0*xiVir[p][sigma]*eta[lambda][nu];
2456+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
2457+ temp = 4.0*xiVir[p][sigma]*eta[mu][nu]
2458+ -1.0*xiVir[p][mu]*eta[nu][sigma]
2459+ -1.0*xiVir[p][nu]*eta[mu][sigma];
2460+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
2461+ }
2462+
2463+ if(lambda!=sigma){
2464+ if(isMoPOcc){
2465+ int p = numberOcc - (moP+1);
2466+ temp = 4.0*xiOcc[p][nu]*eta[sigma][lambda]
2467+ -1.0*xiOcc[p][sigma]*eta[nu][lambda]
2468+ -1.0*xiOcc[p][lambda]*eta[nu][sigma];
2469+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
2470+ temp = 4.0*xiOcc[p][lambda]*eta[mu][nu]
2471+ -1.0*xiOcc[p][mu]*eta[lambda][nu]
2472+ -1.0*xiOcc[p][nu]*eta[lambda][mu];
2473+ value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
2474+ }
2475+ else{
2476+ int p = moP - numberOcc;
2477+ temp = 4.0*xiVir[p][nu]*eta[sigma][lambda]
2478+ -1.0*xiVir[p][sigma]*eta[lambda][nu]
2479+ -1.0*xiVir[p][lambda]*eta[sigma][nu];
2480+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
2481+ temp = 4.0*xiVir[p][lambda]*eta[mu][nu]
2482+ -1.0*xiVir[p][mu]*eta[nu][lambda]
2483+ -1.0*xiVir[p][nu]*eta[mu][lambda];
2484+ value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
2485+ }
2486+ }
2487+
2488+ if(mu!=nu){
2489+ if(isMoPOcc){
2490+ int p = numberOcc - (moP+1);
2491+ temp = 4.0*xiOcc[p][mu]*eta[lambda][sigma]
2492+ -1.0*xiOcc[p][lambda]*eta[mu][sigma]
2493+ -1.0*xiOcc[p][sigma]*eta[mu][lambda];
2494+ value += twoElecInt*this->fockMatrix[moI][nu]*temp;
2495+ temp = 4.0*xiOcc[p][sigma]*eta[nu][mu]
2496+ -1.0*xiOcc[p][nu]*eta[sigma][mu]
2497+ -1.0*xiOcc[p][mu]*eta[sigma][nu];
2498+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
2499+ }
2500+ else{
2501+ int p = moP - numberOcc;
2502+ temp = 4.0*xiVir[p][mu]*eta[lambda][sigma]
2503+ -1.0*xiVir[p][lambda]*eta[sigma][mu]
2504+ -1.0*xiVir[p][sigma]*eta[lambda][mu];
2505+ value += twoElecInt*this->fockMatrix[moI][nu]*temp;
2506+ temp = 4.0*xiVir[p][sigma]*eta[nu][mu]
2507+ -1.0*xiVir[p][nu]*eta[mu][sigma]
2508+ -1.0*xiVir[p][mu]*eta[nu][sigma];
2509+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
2510+ }
2511+ }
2512+
2513+ if(mu!=nu && lambda!=sigma){
2514+ if(isMoPOcc){
2515+ int p = numberOcc - (moP+1);
2516+ temp = 4.0*xiOcc[p][mu]*eta[sigma][lambda]
2517+ -1.0*xiOcc[p][sigma]*eta[mu][lambda]
2518+ -1.0*xiOcc[p][lambda]*eta[mu][sigma];
2519+ value += twoElecInt*this->fockMatrix[moI][nu]*temp;
2520+ temp = 4.0*xiOcc[p][lambda]*eta[nu][mu]
2521+ -1.0*xiOcc[p][nu]*eta[lambda][mu]
2522+ -1.0*xiOcc[p][mu]*eta[lambda][nu];
2523+ value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
2524+ }
2525+ else{
2526+ int p = moP - numberOcc;
2527+ temp = 4.0*xiVir[p][mu]*eta[sigma][lambda]
2528+ -1.0*xiVir[p][sigma]*eta[lambda][mu]
2529+ -1.0*xiVir[p][lambda]*eta[sigma][mu];
2530+ value += twoElecInt*this->fockMatrix[moI][nu]*temp;
2531+ temp = 4.0*xiVir[p][lambda]*eta[nu][mu]
2532+ -1.0*xiVir[p][nu]*eta[mu][lambda]
2533+ -1.0*xiVir[p][mu]*eta[nu][lambda];
2534+ value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
2535+ }
2536+ }
2537+
2538+ }
2539+ }
2540+ }
2541+ }
2542+ }
2543+ else{
2544+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
2545+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
2546+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
2547+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
2548+ double twoElecInt = 0.0;
2549+ if(mu==nu && lambda==sigma){
2550+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2551+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
2552+ twoElecInt = this->GetCoulombInt(orbitalMu,
2553+ orbitalLambda,
2554+ atomA);
2555+ }
2556+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
2557+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
2558+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
2559+ twoElecInt = this->GetExchangeInt(orbitalMu,
2560+ orbitalNu,
2561+ atomA);
2562+ }
2563+ else{
2564+ twoElecInt = 0.0;
2565+ }
2566+
2567+ double temp = 0.0;
2568+ if(isMoPOcc){
2569+ int p = numberOcc - (moP+1);
2570+ temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
2571+ -1.0*xiOcc[p][lambda]*eta[nu][sigma]
2572+ -1.0*xiOcc[p][sigma]*eta[nu][lambda];
2573+ }
2574+ else{
2575+ int p = moP - numberOcc;
2576+ temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
2577+ -1.0*xiVir[p][lambda]*eta[sigma][nu]
2578+ -1.0*xiVir[p][sigma]*eta[lambda][nu];
2579+ }
2580+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
2581+ }
2582+ }
2583+ }
2584+ }
2585+ }
2586+ }
2587+ }
2588+ return value;
2589+}
2590+
24372591 // see common term in eqs. (45) and (46) in [PT_1996],
24382592 // that is, 4.0(ij|kl) - (ik|jl) - (il|jk).
24392593 double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
@@ -6512,4 +6666,3 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
65126666 }
65136667
65146668 }
6515-
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -101,6 +101,11 @@ protected:
101101 double const* const* fockMatrix,
102102 double const* const* gammaAB) const;
103103 virtual void CalcCISMatrix(double** matrixCIS) const;
104+ virtual double GetSmallQElement(int moI,
105+ int moP,
106+ double const* const* xiOcc,
107+ double const* const* xiVir,
108+ double const* const* eta) const;
104109 virtual double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
105110 private:
106111 std::string errorMessageMultipoleA;
@@ -376,11 +381,6 @@ private:
376381 int indexAtomA,
377382 int indexAtomB,
378383 double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
379- void CalcForceExcitedOverlapAOsPart(double* force,
380- int elecStateIndex,
381- int indexAtomA,
382- int indexAtomB,
383- double const* const* const* diatomicOverlapAOs1stDerivs) const;
384384 void CalcForceExcitedTwoElecPart(double* force,
385385 int elecStateIndex,
386386 int indexAtomA,
--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -140,8 +140,6 @@ void ZindoS::SetMessages(){
140140 this->errorMessageDavidsonNotConverged = "Error in zindo::ZindoS::DoCISDavidson: Davidson did not met convergence criterion. \n";
141141 this->errorMessageDavidsonMaxIter = "Davidson loop reaches max_iter=";
142142 this->errorMessageDavidsonMaxDim = "Dimension of the expansion vectors reaches max_dim=";
143- this->errorMessageCalcForceNotGroundState
144- = "Error in zindo::ZindoS::CalcForce: Only ground state is enable in ZindoS.";
145143 this->errorMessageElecState = "Electronic State = ";
146144 this->errorMessageGetElectronicEnergyEnergyNotCalculated
147145 = "Error in zindo::ZindoS::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
@@ -2704,48 +2702,46 @@ void ZindoS::CalcZMatrixForce(const vector<int>& elecStates){
27042702 this->CalcKRDagerGammaRInvMatrix(kRDagerGammaRInv, nonRedundantQIndeces,redundantQIndeces);
27052703 int groundState=0;
27062704 for(int n=0; n<elecStates.size(); n++){
2707- if(groundState < elecStates[n]){
2708- int exciteState = elecStates[n]-1;
2709- this->CalcDeltaVector(delta, exciteState);
2710- this->CalcXiMatrices(xiOcc, xiVir, exciteState, transposedFockMatrix);
2711- this->CalcQVector(q,
2712- delta,
2713- xiOcc,
2714- xiVir,
2715- this->etaMatrixForce[n],
2716- nonRedundantQIndeces,
2717- redundantQIndeces);
2718- this->CalcAuxiliaryVector(y, q, kRDagerGammaRInv, nonRedundantQIndeces, redundantQIndeces);
2719- // solve (54) in [PT_1996]
2720- MolDS_wrappers::Lapack::GetInstance()->Dsysv(gammaNRMinusKNR,
2721- y,
2722- nonRedundantQIndeces.size());
2723- // calculate each element of Z matrix.
2724- stringstream ompErrors;
2705+ if(elecStates[n] <= groundState){continue;}
2706+ int exciteState = elecStates[n]-1;
2707+ this->CalcDeltaVector(delta, exciteState);
2708+ this->CalcXiMatrices(xiOcc, xiVir, exciteState, transposedFockMatrix);
2709+ this->CalcQVector(q,
2710+ delta,
2711+ xiOcc,
2712+ xiVir,
2713+ this->etaMatrixForce[n],
2714+ nonRedundantQIndeces,
2715+ redundantQIndeces);
2716+ this->CalcAuxiliaryVector(y, q, kRDagerGammaRInv, nonRedundantQIndeces, redundantQIndeces);
2717+ // solve (54) in [PT_1996]
2718+ MolDS_wrappers::Lapack::GetInstance()->Dsysv(gammaNRMinusKNR,
2719+ y,
2720+ nonRedundantQIndeces.size());
2721+ // calculate each element of Z matrix.
2722+ stringstream ompErrors;
27252723 #pragma omp parallel for schedule(auto)
2726- for(int mu=0; mu<this->molecule->GetTotalNumberAOs(); mu++){
2727- try{
2728- for(int nu=0; nu<this->molecule->GetTotalNumberAOs(); nu++){
2729- this->zMatrixForce[n][mu][nu] = this->GetZMatrixForceElement(
2730- y,
2731- q,
2732- transposedFockMatrix,
2733- nonRedundantQIndeces,
2734- redundantQIndeces,
2735- mu,
2736- nu);
2737- }
2738- }
2739- catch(MolDSException ex){
2740-#pragma omp critical
2741- ompErrors << ex.what() << endl ;
2724+ for(int mu=0; mu<this->molecule->GetTotalNumberAOs(); mu++){
2725+ try{
2726+ for(int nu=0; nu<this->molecule->GetTotalNumberAOs(); nu++){
2727+ this->zMatrixForce[n][mu][nu] = this->GetZMatrixForceElement(
2728+ y,
2729+ q,
2730+ transposedFockMatrix,
2731+ nonRedundantQIndeces,
2732+ redundantQIndeces,
2733+ mu,
2734+ nu);
27422735 }
27432736 }
2744- // Exception throwing for omp-region
2745- if(!ompErrors.str().empty()){
2746- throw MolDSException(ompErrors.str());
2737+ catch(MolDSException ex){
2738+#pragma omp critical
2739+ ompErrors << ex.what() << endl ;
27472740 }
2748-
2741+ }
2742+ // Exception throwing for omp-region
2743+ if(!ompErrors.str().empty()){
2744+ throw MolDSException(ompErrors.str());
27492745 }
27502746 }
27512747 }
@@ -3005,7 +3001,7 @@ void ZindoS::CalcQVector(double* q,
30053001 */
30063002 }
30073003
3008-// see (18) in [PT_1977]
3004+// see (18) in [PT_1997]
30093005 double ZindoS::GetSmallQElement(int moI,
30103006 int moP,
30113007 double const* const* xiOcc,
@@ -3026,117 +3022,43 @@ double ZindoS::GetSmallQElement(int moI,
30263022 int lastAOIndexB = atomB.GetLastAOIndex();
30273023
30283024 if(A!=B){
3025+ double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
30293026 for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3030- for(int nu=mu; nu<=lastAOIndexA; nu++){
3031- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3032- for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
3033- double twoElecInt = 0.0;
3034- twoElecInt = this->twoElecTwoCore[A]
3035- [B]
3036- [mu-firstAOIndexA]
3037- [nu-firstAOIndexA]
3038- [lambda-firstAOIndexB]
3039- [sigma-firstAOIndexB];
3040- double temp = 0.0;
3041- if(isMoPOcc){
3042- int p = numberOcc - (moP+1);
3043- temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
3044- -1.0*xiOcc[p][lambda]*eta[nu][sigma]
3045- -1.0*xiOcc[p][sigma]*eta[nu][lambda];
3046- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3047- temp = 4.0*xiOcc[p][sigma]*eta[mu][nu]
3048- -1.0*xiOcc[p][mu]*eta[sigma][nu]
3049- -1.0*xiOcc[p][nu]*eta[sigma][mu];
3050- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3051- }
3052- else{
3053- int p = moP - numberOcc;
3054- temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
3055- -1.0*xiVir[p][lambda]*eta[sigma][nu]
3056- -1.0*xiVir[p][sigma]*eta[lambda][nu];
3057- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3058- temp = 4.0*xiVir[p][sigma]*eta[mu][nu]
3059- -1.0*xiVir[p][mu]*eta[nu][sigma]
3060- -1.0*xiVir[p][nu]*eta[mu][sigma];
3061- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3062- }
3063-
3064- if(lambda!=sigma){
3065- if(isMoPOcc){
3066- int p = numberOcc - (moP+1);
3067- temp = 4.0*xiOcc[p][nu]*eta[sigma][lambda]
3068- -1.0*xiOcc[p][sigma]*eta[nu][lambda]
3069- -1.0*xiOcc[p][lambda]*eta[nu][sigma];
3070- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3071- temp = 4.0*xiOcc[p][lambda]*eta[mu][nu]
3072- -1.0*xiOcc[p][mu]*eta[lambda][nu]
3073- -1.0*xiOcc[p][nu]*eta[lambda][mu];
3074- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
3075- }
3076- else{
3077- int p = moP - numberOcc;
3078- temp = 4.0*xiVir[p][nu]*eta[sigma][lambda]
3079- -1.0*xiVir[p][sigma]*eta[lambda][nu]
3080- -1.0*xiVir[p][lambda]*eta[sigma][nu];
3081- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3082- temp = 4.0*xiVir[p][lambda]*eta[mu][nu]
3083- -1.0*xiVir[p][mu]*eta[nu][lambda]
3084- -1.0*xiVir[p][nu]*eta[mu][lambda];
3085- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
3086- }
3087- }
3088-
3089- if(mu!=nu){
3090- if(isMoPOcc){
3091- int p = numberOcc - (moP+1);
3092- temp = 4.0*xiOcc[p][mu]*eta[lambda][sigma]
3093- -1.0*xiOcc[p][lambda]*eta[mu][sigma]
3094- -1.0*xiOcc[p][sigma]*eta[mu][lambda];
3095- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
3096- temp = 4.0*xiOcc[p][sigma]*eta[nu][mu]
3097- -1.0*xiOcc[p][nu]*eta[sigma][mu]
3098- -1.0*xiOcc[p][mu]*eta[sigma][nu];
3099- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3100- }
3101- else{
3102- int p = moP - numberOcc;
3103- temp = 4.0*xiVir[p][mu]*eta[lambda][sigma]
3104- -1.0*xiVir[p][lambda]*eta[sigma][mu]
3105- -1.0*xiVir[p][sigma]*eta[lambda][mu];
3106- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
3107- temp = 4.0*xiVir[p][sigma]*eta[nu][mu]
3108- -1.0*xiVir[p][nu]*eta[mu][sigma]
3109- -1.0*xiVir[p][mu]*eta[nu][sigma];
3110- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3111- }
3112- }
3113-
3114- if(mu!=nu && lambda!=sigma){
3115- if(isMoPOcc){
3116- int p = numberOcc - (moP+1);
3117- temp = 4.0*xiOcc[p][mu]*eta[sigma][lambda]
3118- -1.0*xiOcc[p][sigma]*eta[mu][lambda]
3119- -1.0*xiOcc[p][lambda]*eta[mu][sigma];
3120- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
3121- temp = 4.0*xiOcc[p][lambda]*eta[nu][mu]
3122- -1.0*xiOcc[p][nu]*eta[lambda][mu]
3123- -1.0*xiOcc[p][mu]*eta[lambda][nu];
3124- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
3125- }
3126- else{
3127- int p = moP - numberOcc;
3128- temp = 4.0*xiVir[p][mu]*eta[sigma][lambda]
3129- -1.0*xiVir[p][sigma]*eta[lambda][mu]
3130- -1.0*xiVir[p][lambda]*eta[sigma][mu];
3131- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
3132- temp = 4.0*xiVir[p][lambda]*eta[nu][mu]
3133- -1.0*xiVir[p][nu]*eta[mu][lambda]
3134- -1.0*xiVir[p][mu]*eta[nu][lambda];
3135- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
3136- }
3137- }
3138-
3139- }
3027+ const OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3028+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3029+ const OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3030+ double twoElecInt = 0.0;
3031+ twoElecInt = this->GetNishimotoMatagaTwoEleInt(atomA,
3032+ orbitalMu,
3033+ atomB,
3034+ orbitalLambda,
3035+ rAB);
3036+ double temp = 0.0;
3037+ if(isMoPOcc){
3038+ int p = numberOcc - (moP+1);
3039+ temp = 4.0*xiOcc[p][mu] *eta[lambda][lambda]
3040+ -1.0*xiOcc[p][lambda]*eta[mu][lambda]
3041+ -1.0*xiOcc[p][lambda]*eta[mu][lambda];
3042+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3043+
3044+ temp = 4.0*xiOcc[p][lambda]*eta[mu][mu]
3045+ -1.0*xiOcc[p][mu] *eta[lambda][mu]
3046+ -1.0*xiOcc[p][mu] *eta[lambda][mu];
3047+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3048+
3049+ }
3050+ else{
3051+ int p = moP - numberOcc;
3052+ temp = 4.0*xiVir[p][mu] *eta[lambda][lambda]
3053+ -1.0*xiVir[p][lambda]*eta[lambda][mu]
3054+ -1.0*xiVir[p][lambda]*eta[lambda][mu];
3055+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3056+
3057+ temp = 4.0*xiVir[p][lambda]*eta[mu][mu]
3058+ -1.0*xiVir[p][mu] *eta[mu][lambda]
3059+ -1.0*xiVir[p][mu] *eta[mu][lambda];
3060+ value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
3061+
31403062 }
31413063 }
31423064 }
@@ -3845,16 +3767,11 @@ void ZindoS::CalcDiatomicTwoElecTwoCore1stDerivatives(double*** matrix,
38453767 // elecStates is indeces of the electroinc eigen states.
38463768 // The index = 0 means electronic ground state.
38473769 void ZindoS::CalcForce(const vector<int>& elecStates){
3848- int elecState = elecStates[0];
3849- int groundState = 0;
3850- if(elecState != groundState){
3851- stringstream ss;
3852- ss << this->errorMessageCalcForceNotGroundState;
3853- ss << this->errorMessageElecState << elecState << "\n";
3854- throw MolDSException(ss.str());
3855- }
3856-
38573770 this->CheckMatrixForce(elecStates);
3771+ if(this->RequiresExcitedStatesForce(elecStates)){
3772+ this->CalcEtaMatrixForce(elecStates);
3773+ this->CalcZMatrixForce(elecStates);
3774+ }
38583775 stringstream ompErrors;
38593776 #pragma omp parallel
38603777 {
@@ -3962,6 +3879,13 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
39623879 a,
39633880 b,
39643881 diatomicTwoElecTwoCore1stDerivs);
3882+ // two electron part (excited states)
3883+ double forceExcitedTwoElecPart[CartesianType_end] = {0.0,0.0,0.0};
3884+ this->CalcForceExcitedTwoElecPart(forceExcitedTwoElecPart,
3885+ n,
3886+ a,
3887+ b,
3888+ diatomicTwoElecTwoCore1stDerivs);
39653889 }
39663890 } // end of for(int b)
39673891 } // end of for(int a)
@@ -4112,6 +4036,66 @@ void ZindoS::CalcForceExcitedElecCoreAttractionPart(double* force,
41124036 }
41134037 }
41144038
4039+void ZindoS::CalcForceExcitedOverlapAOsPart(double* force,
4040+ int elecStateIndex,
4041+ int indexAtomA,
4042+ int indexAtomB,
4043+ double const* const* const* diatomicOverlapAOs1stDerivs) const{
4044+ const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
4045+ const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
4046+ int firstAOIndexA = atomA.GetFirstAOIndex();
4047+ int firstAOIndexB = atomB.GetFirstAOIndex();
4048+ int lastAOIndexA = atomA.GetLastAOIndex();
4049+ int lastAOIndexB = atomB.GetLastAOIndex();
4050+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
4051+ for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
4052+ double bondParameter = atomA.GetBondingParameter(
4053+ this->theory,
4054+ atomA.GetValence(mu-firstAOIndexA))
4055+ +atomB.GetBondingParameter(
4056+ this->theory,
4057+ atomB.GetValence(nu-firstAOIndexB));
4058+ bondParameter *= 0.5;
4059+ for(int i=0; i<CartesianType_end; i++){
4060+ force[i] += -1.0
4061+ *this->zMatrixForce[elecStateIndex][mu][nu]
4062+ *bondParameter
4063+ *diatomicOverlapAOs1stDerivs[mu-firstAOIndexA][nu-firstAOIndexB][i];
4064+ }
4065+ }
4066+ }
4067+}
4068+
4069+void ZindoS::CalcForceExcitedTwoElecPart(double* force,
4070+ int elecStateIndex,
4071+ int indexAtomA,
4072+ int indexAtomB,
4073+ double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
4074+ const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
4075+ const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
4076+ int firstAOIndexA = atomA.GetFirstAOIndex();
4077+ int firstAOIndexB = atomB.GetFirstAOIndex();
4078+ int lastAOIndexA = atomA.GetLastAOIndex();
4079+ int lastAOIndexB = atomB.GetLastAOIndex();
4080+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
4081+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
4082+ for(int i=0; i<CartesianType_end; i++){
4083+ force[i] -= this->zMatrixForce[elecStateIndex][mu][mu]
4084+ *this->orbitalElectronPopulation[lambda][lambda]
4085+ *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
4086+ [lambda-firstAOIndexB]
4087+ [i];
4088+ force[i] += 0.50
4089+ *this->zMatrixForce[elecStateIndex][mu][lambda]
4090+ *this->orbitalElectronPopulation[mu][lambda]
4091+ *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
4092+ [lambda-firstAOIndexB]
4093+ [i];
4094+ }
4095+ }
4096+ }
4097+}
4098+
41154099
41164100 }
41174101
--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -95,11 +95,11 @@ protected:
9595 double const* const* eta,
9696 const std::vector<MoIndexPair>& nonRedundantQIndeces,
9797 const std::vector<MoIndexPair>& redundantQIndeces) const;
98- double GetSmallQElement(int moI,
99- int moP,
100- double const* const* xiOcc,
101- double const* const* xiVir,
102- double const* const* eta) const;
98+ virtual double GetSmallQElement(int moI,
99+ int moP,
100+ double const* const* xiOcc,
101+ double const* const* xiVir,
102+ double const* const* eta) const;
103103 void CalcXiMatrices(double** xiOcc,
104104 double** xiVir,
105105 int exciteState,
@@ -122,6 +122,11 @@ protected:
122122 void CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
123123 const std::vector<MoIndexPair>& nonRedundantQIndeces,
124124 const std::vector<MoIndexPair>& redundantQIndeces) const;
125+ void CalcForceExcitedOverlapAOsPart(double* force,
126+ int elecStateIndex,
127+ int indexAtomA,
128+ int indexAtomB,
129+ double const* const* const* diatomicOverlapAOs1stDerivs) const;
125130 /*** end from MNDO ***/
126131 virtual void SetMessages();
127132 virtual void SetEnableAtomTypes();
@@ -199,7 +204,6 @@ protected:
199204 int GetActiveVirIndex(const MolDS_base::Molecule& molecule, int matrixCISIndex) const;
200205 void CheckMatrixForce(const std::vector<int>& elecStates);
201206 private:
202- std::string errorMessageCalcForceNotGroundState;
203207 std::string errorMessageElecState;
204208 std::string errorMessageNishimotoMataga;
205209 std::string errorMessageDavidsonMaxIter;
@@ -247,6 +251,11 @@ private:
247251 int indexAtomA,
248252 int indexAtomB,
249253 double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
254+ void CalcForceExcitedTwoElecPart(double* force,
255+ int elecStateIndex,
256+ int indexAtomA,
257+ int indexAtomB,
258+ double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
250259 void CalcFreeExcitonEnergies(double** freeExcitonEnergiesCIS,
251260 const MolDS_base::Molecule& molecule,
252261 double const* energiesMO,