OSDN -- 開發和下載開源軟體
繁體中文翻譯狀態 translation status for zh-tw
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

提交

標籤
無標籤

Frequently used words (click to add to your profile)

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

Commit MetaInfo

修訂16c75ec423ee98015e09a9dc6f6652e13eaf3aeb (tree)
時間2013-04-25 18:02:36
作者Mikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Some methods for excited forces in Mndo are moved to ZindoS temporary. #31221

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

Change Summary

差異

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -1505,6 +1505,14 @@ double Cndo2::GetFockOffDiagElement(const Atom& atomA,
15051505 return value;
15061506 }
15071507
1508+void Cndo2::TransposeFockMatrixMatrix(double** transposedFockMatrix) const{
1509+ for(int i=0; i<this->molecule->GetTotalNumberAOs(); i++){
1510+ for(int j=0; j<this->molecule->GetTotalNumberAOs(); j++){
1511+ transposedFockMatrix[j][i] = this->fockMatrix[i][j];
1512+ }
1513+ }
1514+}
1515+
15081516 void Cndo2::CalcOrbitalElectronPopulation(double** orbitalElectronPopulation,
15091517 const Molecule& molecule,
15101518 double const* const* fockMatrix) const{
--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -184,6 +184,7 @@ protected:
184184 double const* const* orbitalElectronPopulation,
185185 double const* const* const* const* const* const* twoElecTwoCore,
186186 bool isGuess) const;
187+ void TransposeFockMatrixMatrix(double** transposedFockMatrix) const;
187188 virtual void CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs,
188189 const MolDS_base_atoms::Atom& atomA,
189190 const MolDS_base_atoms::Atom& atomB) const;
--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -62,10 +62,6 @@ Mndo::Mndo() : MolDS_zindo::ZindoS(){
6262 this->SetEnableAtomTypes();
6363 // private variables
6464 this->heatsFormation = 0.0;
65- this->zMatrixForceElecStatesNum = 0;
66- this->etaMatrixForceElecStatesNum = 0;
67- this->zMatrixForce = NULL;
68- this->etaMatrixForce = NULL;
6965 //this->OutputLog("Mndo created\n");
7066 }
7167
@@ -78,14 +74,6 @@ Mndo::~Mndo(){
7874 dxy,
7975 dxy,
8076 dxy);
81- MallocerFreer::GetInstance()->Free<double>(&this->zMatrixForce,
82- this->zMatrixForceElecStatesNum,
83- this->molecule->GetTotalNumberAOs(),
84- this->molecule->GetTotalNumberAOs());
85- MallocerFreer::GetInstance()->Free<double>(&this->etaMatrixForce,
86- this->etaMatrixForceElecStatesNum,
87- this->molecule->GetTotalNumberAOs(),
88- this->molecule->GetTotalNumberAOs());
8977 MallocerFreer::GetInstance()->Free<double>(&this->normalForceConstants,
9078 CartesianType_end*molecule->GetNumberAtoms());
9179 MallocerFreer::GetInstance()->Free<double>(&this->normalModes,
@@ -906,125 +894,6 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
906894 % this->messageDoneCalcCISMatrix.c_str());
907895 }
908896
909-void Mndo::CheckZMatrixForce(const vector<int>& elecStates){
910- // malloc or initialize Z matrix
911- if(this->zMatrixForce == NULL){
912- MallocerFreer::GetInstance()->Malloc<double>(&this->zMatrixForce,
913- elecStates.size(),
914- this->molecule->GetTotalNumberAOs(),
915- this->molecule->GetTotalNumberAOs());
916- this->zMatrixForceElecStatesNum = elecStates.size();
917- }
918- else{
919- MallocerFreer::GetInstance()->
920- Initialize<double>(this->zMatrixForce,
921- elecStates.size(),
922- this->molecule->GetTotalNumberAOs(),
923- this->molecule->GetTotalNumberAOs());
924- }
925-}
926-
927-void Mndo::CheckEtaMatrixForce(const vector<int>& elecStates){
928- // malloc or initialize eta matrix
929- if(this->etaMatrixForce == NULL){
930- MallocerFreer::GetInstance()->Malloc<double>(&this->etaMatrixForce,
931- elecStates.size(),
932- this->molecule->GetTotalNumberAOs(),
933- this->molecule->GetTotalNumberAOs());
934- this->etaMatrixForceElecStatesNum = elecStates.size();
935- }
936- else{
937- MallocerFreer::GetInstance()->
938- Initialize<double>(this->etaMatrixForce,
939- elecStates.size(),
940- this->molecule->GetTotalNumberAOs(),
941- this->molecule->GetTotalNumberAOs());
942- }
943-}
944-
945-// see variable Q-vector in [PT_1996, PT_1997]
946-void Mndo::CalcActiveSetVariablesQ(vector<MoIndexPair>* nonRedundantQIndeces,
947- vector<MoIndexPair>* redundantQIndeces,
948- int numberActiveOcc,
949- int numberActiveVir) const{
950- int numberAOs = this->molecule->GetTotalNumberAOs();
951- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
952- for(int moI=0; moI<numberOcc; moI++){
953- bool isMoICIMO = numberOcc-numberActiveOcc<=moI ? true : false;
954- for(int moJ=numberOcc; moJ<numberAOs; moJ++){
955- bool isMoJCIMO = moJ<numberOcc+numberActiveVir ? true : false;
956- MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
957- nonRedundantQIndeces->push_back(moIndexPair);
958- }
959- }
960- for(int moI=numberOcc-numberActiveOcc; moI<numberOcc; moI++){
961- bool isMoICIMO = true;
962- for(int moJ=moI; moJ<numberOcc; moJ++){
963- bool isMoJCIMO = true;
964- MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
965- redundantQIndeces->push_back(moIndexPair);
966- }
967- }
968- for(int moI=numberOcc; moI<numberOcc+numberActiveVir; moI++){
969- bool isMoICIMO = true;
970- for(int moJ=moI; moJ<numberOcc+numberActiveVir; moJ++){
971- bool isMoJCIMO = true;
972- MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
973- redundantQIndeces->push_back(moIndexPair);
974- }
975- }
976-}
977-
978-void Mndo::MallocTempMatrixForZMatrix(double** delta,
979- double** q,
980- double*** gammaNRMinusKNR,
981- double*** kRDag,
982- double** y,
983- double*** transposedFockMatrix,
984- double*** xiOcc,
985- double*** xiVir,
986- int sizeQNR,
987- int sizeQR) const{
988- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
989- int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
990- int numberActiveMO = numberActiveOcc + numberActiveVir;
991- int numberAOs = this->molecule->GetTotalNumberAOs();
992- MallocerFreer::GetInstance()->Malloc<double>(delta, numberActiveMO);
993- MallocerFreer::GetInstance()->Malloc<double>(q, sizeQNR+sizeQR);
994- MallocerFreer::GetInstance()->Malloc<double>(gammaNRMinusKNR, sizeQNR, sizeQNR);
995- MallocerFreer::GetInstance()->Malloc<double>(kRDag, sizeQNR, sizeQR);
996- MallocerFreer::GetInstance()->Malloc<double>(y, sizeQNR);
997- MallocerFreer::GetInstance()->Malloc<double>(transposedFockMatrix,
998- numberAOs,
999- numberAOs);
1000- MallocerFreer::GetInstance()->Malloc<double>(xiOcc, numberActiveOcc,numberAOs);
1001- MallocerFreer::GetInstance()->Malloc<double>(xiVir,numberActiveVir,numberAOs);
1002-}
1003-
1004-void Mndo::FreeTempMatrixForZMatrix(double** delta,
1005- double** q,
1006- double*** gammaNRMinusKNR,
1007- double*** kRDag,
1008- double** y,
1009- double*** transposedFockMatrix,
1010- double*** xiOcc,
1011- double*** xiVir,
1012- int sizeQNR,
1013- int sizeQR) const{
1014- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
1015- int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
1016- int numberActiveMO = numberActiveOcc + numberActiveVir;
1017- int numberAOs = this->molecule->GetTotalNumberAOs();
1018- MallocerFreer::GetInstance()->Free<double>(delta, numberActiveMO);
1019- MallocerFreer::GetInstance()->Free<double>(q, sizeQNR+sizeQR);
1020- MallocerFreer::GetInstance()->Free<double>(gammaNRMinusKNR, sizeQNR, sizeQNR);
1021- MallocerFreer::GetInstance()->Free<double>(kRDag, sizeQNR, sizeQR);
1022- MallocerFreer::GetInstance()->Free<double>(y, sizeQNR);
1023- MallocerFreer::GetInstance()->Free<double>(transposedFockMatrix, numberAOs, numberAOs);
1024- MallocerFreer::GetInstance()->Free<double>(xiOcc, numberActiveOcc, numberAOs);
1025- MallocerFreer::GetInstance()->Free<double>(xiVir, numberActiveVir, numberAOs);
1026-}
1027-
1028897 // \epsilon_{r}^{kl} in (1) in [PT_1997].
1029898 // k and l are index of CIS matrix.
1030899 double Mndo::GetCISCoefficientMOEnergy(int k, int l, int r, int numberActiveVir) const{
@@ -1069,974 +938,6 @@ double Mndo::GetCISCoefficientTwoElecIntegral(int k,
1069938 return value;
1070939 }
1071940
1072-// see (40) in [PT_1996]
1073-double Mndo::GetGammaNRElement(int moI, int moJ, int moK, int moL) const{
1074- double value=0.0;
1075- if(moI==moK && moJ==moL){
1076- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1077- double nI = moI<numberOcc ? 2.0 : 0.0;
1078- double nJ = moJ<numberOcc ? 2.0 : 0.0;
1079- value = (this->energiesMO[moJ]-this->energiesMO[moI])/(nJ-nI);
1080- }
1081- return value;
1082-}
1083-
1084-// see (41) & (42) in [PT_1996]
1085-double Mndo::GetGammaRElement(int moI, int moJ, int moK, int moL) const{
1086- double value=0.0;
1087- if(moI==moK && moJ==moL){
1088- value = moI==moJ ? 1.0 : this->energiesMO[moJ]-this->energiesMO[moI];
1089- }
1090- return value;
1091-}
1092-
1093-// see (43) in [PT_1996]
1094-double Mndo::GetNNRElement(int moI, int moJ, int moK, int moL) const{
1095- double value=0.0;
1096- if(moI==moK && moJ==moL){
1097- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1098- double nI = moI<numberOcc ? 2.0 : 0.0;
1099- double nJ = moJ<numberOcc ? 2.0 : 0.0;
1100- value = (nJ-nI);
1101- }
1102- return value;
1103-}
1104-
1105-// see (44) in [PT_1996]
1106-double Mndo::GetNRElement(int moI, int moJ, int moK, int moL) const{
1107- double value=0.0;
1108- if(moI==moK && moJ==moL){
1109- value = 1.0;
1110- }
1111- return value;
1112-}
1113-
1114-// see (44) in [PT_1996]
1115-double Mndo::GetKNRElement(int moI, int moJ, int moK, int moL) const{
1116- double value=0.0;
1117- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1118- int nI = moI<numberOcc ? 2 : 0;
1119- int nJ = moJ<numberOcc ? 2 : 0;
1120- int nK = moK<numberOcc ? 2 : 0;
1121- int nL = moL<numberOcc ? 2 : 0;
1122-
1123- if(nI!=nJ && nK!=nL){
1124- value = this->GetAuxiliaryKNRKRElement(moI, moJ, moK, moL);
1125- }
1126- //See (24) in [DL_1990] about "0.5" multiplied to "GetKNRElement".
1127- return 0.5*value;
1128-}
1129-
1130-// Dager of (45) in [PT_1996]. Note taht the (45) is real number.
1131-double Mndo::GetKRDagerElement(int moI, int moJ, int moK, int moL) const{
1132- return this->GetKRElement(moK, moL, moI, moJ);
1133-}
1134-
1135-// see (45) in [PT_1996]
1136-double Mndo::GetKRElement(int moI, int moJ, int moK, int moL) const{
1137- double value=0.0;
1138- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1139- int nI = moI<numberOcc ? 2 : 0;
1140- int nJ = moJ<numberOcc ? 2 : 0;
1141- int nK = moK<numberOcc ? 2 : 0;
1142- int nL = moL<numberOcc ? 2 : 0;
1143-
1144- if(nI==nJ && nK!=nL){
1145- value = this->GetAuxiliaryKNRKRElement(moI, moJ, moK, moL);
1146- }
1147- //See (24) in [DL_1990] about "0.5" multiplied to "GetKRElement".
1148- return 0.5*value;
1149-}
1150-
1151-double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
1152- double value = 0.0;
1153-
1154- // Fast algorith, but this is not easy to read.
1155- // Slow algorithm is alos written below.
1156- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
1157- const Atom& atomA = *this->molecule->GetAtom(A);
1158- int firstAOIndexA = atomA.GetFirstAOIndex();
1159- int lastAOIndexA = atomA.GetLastAOIndex();
1160-
1161- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1162- int muOffSet = mu - firstAOIndexA;
1163- for(int nu=mu; nu<=lastAOIndexA; nu++){
1164- int nuOffSet = nu - firstAOIndexA;
1165- double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
1166- tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0,
1167- tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0,
1168- tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0;
1169- tmpMN01 = 4.0
1170- *this->fockMatrix[moI][mu]
1171- *this->fockMatrix[moJ][nu];
1172- tmpMN02 = 4.0
1173- *this->fockMatrix[moK][mu]
1174- *this->fockMatrix[moL][nu];
1175- tmpMN03 = this->fockMatrix[moI][mu]
1176- *this->fockMatrix[moK][nu];
1177- tmpMN04 = this->fockMatrix[moJ][mu]
1178- *this->fockMatrix[moL][nu];
1179- tmpMN05 = this->fockMatrix[moI][mu]
1180- *this->fockMatrix[moL][nu];
1181- tmpMN06 = this->fockMatrix[moJ][mu]
1182- *this->fockMatrix[moK][nu];
1183- if(mu != nu){
1184- tmpMN13 = 4.0
1185- *this->fockMatrix[moI][nu]
1186- *this->fockMatrix[moJ][mu];
1187- tmpMN14 = 4.0
1188- *this->fockMatrix[moK][nu]
1189- *this->fockMatrix[moL][mu];
1190- tmpMN15 = this->fockMatrix[moI][nu]
1191- *this->fockMatrix[moK][mu];
1192- tmpMN16 = this->fockMatrix[moJ][nu]
1193- *this->fockMatrix[moL][mu];
1194- tmpMN17 = this->fockMatrix[moI][nu]
1195- *this->fockMatrix[moL][mu];
1196- tmpMN18 = this->fockMatrix[moJ][nu]
1197- *this->fockMatrix[moK][mu];
1198- }
1199-
1200- for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
1201- const Atom& atomB = *this->molecule->GetAtom(B);
1202- int firstAOIndexB = atomB.GetFirstAOIndex();
1203- int lastAOIndexB = atomB.GetLastAOIndex();
1204-
1205- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1206- int lambdaOffSet = lambda - firstAOIndexB;
1207- double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0,
1208- tmpMNL05 = 0.0, tmpMNL06 = 0.0, tmpMNL07 = 0.0, tmpMNL08 = 0.0,
1209- tmpMNL09 = 0.0, tmpMNL10 = 0.0, tmpMNL11 = 0.0, tmpMNL12 = 0.0,
1210- tmpMNL13 = 0.0, tmpMNL14 = 0.0, tmpMNL15 = 0.0, tmpMNL16 = 0.0,
1211- tmpMNL17 = 0.0, tmpMNL18 = 0.0, tmpMNL19 = 0.0, tmpMNL20 = 0.0,
1212- tmpMNL21 = 0.0, tmpMNL22 = 0.0, tmpMNL23 = 0.0, tmpMNL24 = 0.0;
1213- tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda];
1214- tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda];
1215- tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda];
1216- tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda];
1217- tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda];
1218- tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda];
1219- tmpMNL07 = tmpMN01*this->fockMatrix[moL][lambda];
1220- tmpMNL08 = tmpMN02*this->fockMatrix[moJ][lambda];
1221- tmpMNL09 = tmpMN03*this->fockMatrix[moL][lambda];
1222- tmpMNL10 = tmpMN04*this->fockMatrix[moK][lambda];
1223- tmpMNL11 = tmpMN05*this->fockMatrix[moK][lambda];
1224- tmpMNL12 = tmpMN06*this->fockMatrix[moL][lambda];
1225- tmpMNL01 -= tmpMNL03 + tmpMNL06;
1226- tmpMNL04 += tmpMNL05;
1227- tmpMNL08 -= tmpMNL10 + tmpMNL12;
1228- tmpMNL09 += tmpMNL11;
1229- if(mu != nu){
1230- tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda];
1231- tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda];
1232- tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda];
1233- tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda];
1234- tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda];
1235- tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda];
1236- tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda];
1237- tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda];
1238- tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda];
1239- tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda];
1240- tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda];
1241- tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda];
1242- tmpMNL13 -= tmpMNL15 + tmpMNL18;
1243- tmpMNL16 += tmpMNL17;
1244- tmpMNL20 -= tmpMNL22 + tmpMNL24;
1245- tmpMNL21 += tmpMNL23;
1246- tmpMNL01 += tmpMNL13;
1247- tmpMNL02 += tmpMNL14;
1248- tmpMNL04 += tmpMNL16;
1249- tmpMNL07 += tmpMNL19;
1250- tmpMNL08 += tmpMNL20;
1251- tmpMNL09 += tmpMNL21;
1252- }
1253- for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
1254- int sigmaOffSet = sigma - firstAOIndexB;
1255- double tmpValue = 0.0;
1256- tmpValue += tmpMNL01*this->fockMatrix[moL][sigma];
1257- tmpValue += tmpMNL02*this->fockMatrix[moJ][sigma];
1258- tmpValue -= tmpMNL04*this->fockMatrix[moK][sigma];
1259- if(lambda != sigma){
1260- tmpValue += tmpMNL07*this->fockMatrix[moK][sigma];
1261- tmpValue += tmpMNL08*this->fockMatrix[moI][sigma];
1262- tmpValue -= tmpMNL09*this->fockMatrix[moJ][sigma];
1263- }
1264- double gamma = 0.0;
1265- if(A!=B){
1266- gamma = this->twoElecTwoCore[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
1267- }
1268- else{
1269- if(mu==nu && lambda==sigma){
1270- OrbitalType orbitalMu = atomA.GetValence(muOffSet);
1271- OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet);
1272- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
1273- }
1274- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
1275- OrbitalType orbitalMu = atomA.GetValence(muOffSet);
1276- OrbitalType orbitalNu = atomA.GetValence(nuOffSet);
1277- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
1278- }
1279- else{
1280- gamma = 0.0;
1281- }
1282- gamma *= 0.5;
1283- }
1284- value += tmpValue*gamma;
1285- }
1286- }
1287- }
1288- }
1289- }
1290- }
1291- // End of the fast algorith.
1292-
1293- /*
1294- // Algorithm using blas
1295- double** twoElec = NULL;
1296- double* twiceMoIJ = NULL;
1297- double* twiceMoIK = NULL;
1298- double* twiceMoIL = NULL;
1299- double* twiceMoKL = NULL;
1300- double* twiceMoJL = NULL;
1301- double* twiceMoJK = NULL;
1302- double* tmpVector = NULL;
1303- int numAOs = this->molecule->GetTotalNumberAOs();
1304- MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
1305- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
1306- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
1307- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
1308- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
1309- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
1310- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
1311- MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
1312- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
1313- const Atom& atomA = *this->molecule->GetAtom(A);
1314- int firstAOIndexA = atomA.GetFirstAOIndex();
1315- int lastAOIndexA = atomA.GetLastAOIndex();
1316- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1317- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1318- twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
1319- twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
1320- twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
1321- }
1322- }
1323- }
1324-
1325- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
1326- const Atom& atomB = *this->molecule->GetAtom(B);
1327- int firstAOIndexB = atomB.GetFirstAOIndex();
1328- int lastAOIndexB = atomB.GetLastAOIndex();
1329- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1330- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1331- twiceMoKL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
1332- twiceMoJL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
1333- twiceMoJK[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
1334- }
1335- }
1336- }
1337-
1338- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
1339- const Atom& atomA = *this->molecule->GetAtom(A);
1340- int firstAOIndexA = atomA.GetFirstAOIndex();
1341- int lastAOIndexA = atomA.GetLastAOIndex();
1342- for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
1343- const Atom& atomB = *this->molecule->GetAtom(B);
1344- int firstAOIndexB = atomB.GetFirstAOIndex();
1345- int lastAOIndexB = atomB.GetLastAOIndex();
1346- double gamma = 0.0;
1347- if(A!=B){
1348- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1349- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1350- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1351- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1352- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
1353- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
1354- this->twoElecTwoCore[A]
1355- [B]
1356- [mu-firstAOIndexA]
1357- [nu-firstAOIndexA]
1358- [lambda-firstAOIndexB]
1359- [sigma-firstAOIndexB];
1360- }
1361- }
1362- }
1363- }
1364- }
1365- else{
1366- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1367- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1368- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1369- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1370- if(mu==nu && lambda==sigma){
1371- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1372- OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
1373- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
1374- }
1375- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
1376- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1377- OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
1378- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
1379- }
1380- else{
1381- gamma = 0.0;
1382- }
1383- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
1384- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
1385- }
1386- }
1387- }
1388- }
1389- }
1390- }
1391- }
1392- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
1393- twoElec,
1394- twiceMoKL,
1395- tmpVector);
1396- value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, tmpVector);
1397- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
1398- twoElec,
1399- twiceMoJL,
1400- tmpVector);
1401- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, tmpVector);
1402- MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
1403- twoElec,
1404- twiceMoJK,
1405- tmpVector);
1406- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, tmpVector);
1407- MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
1408- MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
1409- MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
1410- MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
1411- MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
1412- MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
1413- MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
1414- MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
1415- // End of algorithm using blas
1416- */
1417-
1418- /*
1419- // Second algorithm using blas.
1420- // This algorithm uses DGEMM.
1421- double** twoElec = NULL;
1422- double* twiceMoIJ = NULL;
1423- double* twiceMoIK = NULL;
1424- double* twiceMoIL = NULL;
1425- double** twiceMoB = NULL;
1426- double** tmpMatrix = NULL;
1427- int numAOs = this->molecule->GetTotalNumberAOs();
1428- MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
1429- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
1430- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
1431- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
1432- MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
1433- MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
1434- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
1435- const Atom& atomA = *this->molecule->GetAtom(A);
1436- int firstAOIndexA = atomA.GetFirstAOIndex();
1437- int lastAOIndexA = atomA.GetLastAOIndex();
1438- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1439- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1440- twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
1441- twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
1442- twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
1443- }
1444- }
1445- }
1446-
1447- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
1448- const Atom& atomB = *this->molecule->GetAtom(B);
1449- int firstAOIndexB = atomB.GetFirstAOIndex();
1450- int lastAOIndexB = atomB.GetLastAOIndex();
1451- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1452- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1453- twiceMoB[0][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
1454- twiceMoB[1][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
1455- twiceMoB[2][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
1456- }
1457- }
1458- }
1459-
1460- for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
1461- const Atom& atomA = *this->molecule->GetAtom(A);
1462- int firstAOIndexA = atomA.GetFirstAOIndex();
1463- int lastAOIndexA = atomA.GetLastAOIndex();
1464- for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
1465- const Atom& atomB = *this->molecule->GetAtom(B);
1466- int firstAOIndexB = atomB.GetFirstAOIndex();
1467- int lastAOIndexB = atomB.GetLastAOIndex();
1468- double gamma = 0.0;
1469- if(A!=B){
1470- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1471- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1472- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1473- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1474- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
1475- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
1476- this->twoElecTwoCore[A]
1477- [B]
1478- [mu-firstAOIndexA]
1479- [nu-firstAOIndexA]
1480- [lambda-firstAOIndexB]
1481- [sigma-firstAOIndexB];
1482- }
1483- }
1484- }
1485- }
1486- }
1487- else{
1488- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1489- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1490- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1491- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1492- if(mu==nu && lambda==sigma){
1493- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1494- OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
1495- gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
1496- }
1497- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
1498- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1499- OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
1500- gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
1501- }
1502- else{
1503- gamma = 0.0;
1504- }
1505- twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
1506- [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
1507- }
1508- }
1509- }
1510- }
1511- }
1512- }
1513- }
1514-
1515- MolDS_wrappers::Blas::GetInstance()->Dgemm(false, true, true,
1516- this->molecule->GetNumberAtoms()*dxy*dxy,
1517- 3,
1518- this->molecule->GetNumberAtoms()*dxy*dxy,
1519- 1.0,
1520- twoElec,
1521- twiceMoB,
1522- 0.0,
1523- tmpMatrix);
1524- value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
1525- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
1526- value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
1527- MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
1528- MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
1529- MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
1530- MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
1531- MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
1532- MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
1533- // End of second algorithm using blas
1534- */
1535-
1536- /*
1537- // slow algorithm
1538- value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL,
1539- *this->molecule,
1540- this->fockMatrix, NULL)
1541- -1.0*this->GetMolecularIntegralElement(moI, moK, moJ, moL,
1542- *this->molecule,
1543- this->fockMatrix, NULL)
1544- -1.0*this->GetMolecularIntegralElement(moI, moL, moJ, moK,
1545- *this->molecule,
1546- this->fockMatrix, NULL);
1547- */
1548- return value;
1549-}
1550-
1551-// see (9) in [PT_1997]
1552-void Mndo::CalcDeltaVector(double* delta, int exciteState) const{
1553- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
1554- int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
1555- int numberActiveMO = numberActiveOcc + numberActiveVir;
1556- MallocerFreer::GetInstance()->Initialize<double>(delta, numberActiveMO);
1557- stringstream ompErrors;
1558-#pragma omp parallel for schedule(auto)
1559- for(int r=0; r<numberActiveMO; r++){
1560- try{
1561- double value = 0.0;
1562- if(r<numberActiveOcc){
1563- // r is active occupied MO
1564- int rr=numberActiveOcc-(r+1);
1565- for(int a=0; a<numberActiveVir; a++){
1566- int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(rr,a);
1567- value -= pow(this->matrixCIS[exciteState][slaterDeterminantIndex],2.0);
1568- }
1569- }
1570- else{
1571- // r is active virtual MO
1572- int rr=r-numberActiveOcc;
1573- for(int i=0; i<numberActiveOcc; i++){
1574- int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,rr);
1575- value += pow(this->matrixCIS[exciteState][slaterDeterminantIndex],2.0);
1576- }
1577- }
1578- delta[r] = value;
1579- }
1580- catch(MolDSException ex){
1581-#pragma omp critical
1582- ompErrors << ex.what() << endl ;
1583- }
1584- }
1585- // Exception throwing for omp-region
1586- if(!ompErrors.str().empty()){
1587- throw MolDSException(ompErrors.str());
1588- }
1589-}
1590-
1591-// see (18) in [PT_1977]
1592-double Mndo::GetSmallQElement(int moI,
1593- int moP,
1594- double const* const* xiOcc,
1595- double const* const* xiVir,
1596- double const* const* eta) const{
1597- double value = 0.0;
1598- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1599- bool isMoPOcc = moP<numberOcc ? true : false;
1600-
1601- for(int A=0; A<molecule->GetNumberAtoms(); A++){
1602- const Atom& atomA = *molecule->GetAtom(A);
1603- int firstAOIndexA = atomA.GetFirstAOIndex();
1604- int lastAOIndexA = atomA.GetLastAOIndex();
1605-
1606- for(int B=A; B<molecule->GetNumberAtoms(); B++){
1607- const Atom& atomB = *molecule->GetAtom(B);
1608- int firstAOIndexB = atomB.GetFirstAOIndex();
1609- int lastAOIndexB = atomB.GetLastAOIndex();
1610-
1611- if(A!=B){
1612- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1613- for(int nu=mu; nu<=lastAOIndexA; nu++){
1614- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1615- for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
1616- double twoElecInt = 0.0;
1617- twoElecInt = this->twoElecTwoCore[A]
1618- [B]
1619- [mu-firstAOIndexA]
1620- [nu-firstAOIndexA]
1621- [lambda-firstAOIndexB]
1622- [sigma-firstAOIndexB];
1623- double temp = 0.0;
1624- if(isMoPOcc){
1625- int p = numberOcc - (moP+1);
1626- temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
1627- -1.0*xiOcc[p][lambda]*eta[nu][sigma]
1628- -1.0*xiOcc[p][sigma]*eta[nu][lambda];
1629- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
1630- temp = 4.0*xiOcc[p][sigma]*eta[mu][nu]
1631- -1.0*xiOcc[p][mu]*eta[sigma][nu]
1632- -1.0*xiOcc[p][nu]*eta[sigma][mu];
1633- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
1634- }
1635- else{
1636- int p = moP - numberOcc;
1637- temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
1638- -1.0*xiVir[p][lambda]*eta[sigma][nu]
1639- -1.0*xiVir[p][sigma]*eta[lambda][nu];
1640- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
1641- temp = 4.0*xiVir[p][sigma]*eta[mu][nu]
1642- -1.0*xiVir[p][mu]*eta[nu][sigma]
1643- -1.0*xiVir[p][nu]*eta[mu][sigma];
1644- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
1645- }
1646-
1647- if(lambda!=sigma){
1648- if(isMoPOcc){
1649- int p = numberOcc - (moP+1);
1650- temp = 4.0*xiOcc[p][nu]*eta[sigma][lambda]
1651- -1.0*xiOcc[p][sigma]*eta[nu][lambda]
1652- -1.0*xiOcc[p][lambda]*eta[nu][sigma];
1653- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
1654- temp = 4.0*xiOcc[p][lambda]*eta[mu][nu]
1655- -1.0*xiOcc[p][mu]*eta[lambda][nu]
1656- -1.0*xiOcc[p][nu]*eta[lambda][mu];
1657- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
1658- }
1659- else{
1660- int p = moP - numberOcc;
1661- temp = 4.0*xiVir[p][nu]*eta[sigma][lambda]
1662- -1.0*xiVir[p][sigma]*eta[lambda][nu]
1663- -1.0*xiVir[p][lambda]*eta[sigma][nu];
1664- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
1665- temp = 4.0*xiVir[p][lambda]*eta[mu][nu]
1666- -1.0*xiVir[p][mu]*eta[nu][lambda]
1667- -1.0*xiVir[p][nu]*eta[mu][lambda];
1668- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
1669- }
1670- }
1671-
1672- if(mu!=nu){
1673- if(isMoPOcc){
1674- int p = numberOcc - (moP+1);
1675- temp = 4.0*xiOcc[p][mu]*eta[lambda][sigma]
1676- -1.0*xiOcc[p][lambda]*eta[mu][sigma]
1677- -1.0*xiOcc[p][sigma]*eta[mu][lambda];
1678- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
1679- temp = 4.0*xiOcc[p][sigma]*eta[nu][mu]
1680- -1.0*xiOcc[p][nu]*eta[sigma][mu]
1681- -1.0*xiOcc[p][mu]*eta[sigma][nu];
1682- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
1683- }
1684- else{
1685- int p = moP - numberOcc;
1686- temp = 4.0*xiVir[p][mu]*eta[lambda][sigma]
1687- -1.0*xiVir[p][lambda]*eta[sigma][mu]
1688- -1.0*xiVir[p][sigma]*eta[lambda][mu];
1689- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
1690- temp = 4.0*xiVir[p][sigma]*eta[nu][mu]
1691- -1.0*xiVir[p][nu]*eta[mu][sigma]
1692- -1.0*xiVir[p][mu]*eta[nu][sigma];
1693- value += twoElecInt*this->fockMatrix[moI][lambda]*temp;
1694- }
1695- }
1696-
1697- if(mu!=nu && lambda!=sigma){
1698- if(isMoPOcc){
1699- int p = numberOcc - (moP+1);
1700- temp = 4.0*xiOcc[p][mu]*eta[sigma][lambda]
1701- -1.0*xiOcc[p][sigma]*eta[mu][lambda]
1702- -1.0*xiOcc[p][lambda]*eta[mu][sigma];
1703- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
1704- temp = 4.0*xiOcc[p][lambda]*eta[nu][mu]
1705- -1.0*xiOcc[p][nu]*eta[lambda][mu]
1706- -1.0*xiOcc[p][mu]*eta[lambda][nu];
1707- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
1708- }
1709- else{
1710- int p = moP - numberOcc;
1711- temp = 4.0*xiVir[p][mu]*eta[sigma][lambda]
1712- -1.0*xiVir[p][sigma]*eta[lambda][mu]
1713- -1.0*xiVir[p][lambda]*eta[sigma][mu];
1714- value += twoElecInt*this->fockMatrix[moI][nu]*temp;
1715- temp = 4.0*xiVir[p][lambda]*eta[nu][mu]
1716- -1.0*xiVir[p][nu]*eta[mu][lambda]
1717- -1.0*xiVir[p][mu]*eta[nu][lambda];
1718- value += twoElecInt*this->fockMatrix[moI][sigma]*temp;
1719- }
1720- }
1721-
1722- }
1723- }
1724- }
1725- }
1726- }
1727- else{
1728- for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
1729- for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
1730- for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
1731- for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
1732- double twoElecInt = 0.0;
1733- if(mu==nu && lambda==sigma){
1734- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1735- OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
1736- twoElecInt = this->GetCoulombInt(orbitalMu,
1737- orbitalLambda,
1738- atomA);
1739- }
1740- else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
1741- OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
1742- OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
1743- twoElecInt = this->GetExchangeInt(orbitalMu,
1744- orbitalNu,
1745- atomA);
1746- }
1747- else{
1748- twoElecInt = 0.0;
1749- }
1750-
1751- double temp = 0.0;
1752- if(isMoPOcc){
1753- int p = numberOcc - (moP+1);
1754- temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
1755- -1.0*xiOcc[p][lambda]*eta[nu][sigma]
1756- -1.0*xiOcc[p][sigma]*eta[nu][lambda];
1757- }
1758- else{
1759- int p = moP - numberOcc;
1760- temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
1761- -1.0*xiVir[p][lambda]*eta[sigma][nu]
1762- -1.0*xiVir[p][sigma]*eta[lambda][nu];
1763- }
1764- value += twoElecInt*this->fockMatrix[moI][mu]*temp;
1765- }
1766- }
1767- }
1768- }
1769- }
1770- }
1771- }
1772- return value;
1773-}
1774-
1775-// see (20) - (23) in [PT_1997]
1776-void Mndo::CalcQVector(double* q,
1777- double const* delta,
1778- double const* const* xiOcc,
1779- double const* const* xiVir,
1780- double const* const* eta,
1781- const vector<MoIndexPair>& nonRedundantQIndeces,
1782- const vector<MoIndexPair>& redundantQIndeces) const{
1783- MallocerFreer::GetInstance()->Initialize<double>(
1784- q,
1785- nonRedundantQIndeces.size()+redundantQIndeces.size());
1786-
1787- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1788- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
1789- stringstream ompErrors;
1790-#pragma omp parallel for schedule(auto)
1791- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1792- try{
1793- int moI = nonRedundantQIndeces[i].moI;
1794- int moJ = nonRedundantQIndeces[i].moJ;
1795- bool isMoICIMO = nonRedundantQIndeces[i].isMoICIMO;
1796- bool isMoJCIMO = nonRedundantQIndeces[i].isMoJCIMO;
1797- if(!isMoICIMO && isMoJCIMO){
1798- q[i] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta);
1799- }
1800- else if(isMoICIMO && !isMoJCIMO){
1801- q[i] = -1.0*this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
1802- }
1803- else if(isMoICIMO && isMoJCIMO){
1804- q[i] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta)
1805- -this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
1806- }
1807- else{
1808- q[i] = 0.0;
1809- }
1810- }
1811- catch(MolDSException ex){
1812-#pragma omp critical
1813- ompErrors << ex.what() << endl ;
1814- }
1815- }
1816- // Exception throwing for omp-region
1817- if(!ompErrors.str().empty()){
1818- throw MolDSException(ompErrors.str());
1819- }
1820-#pragma omp parallel for schedule(auto)
1821- for(int i=0; i<redundantQIndeces.size(); i++){
1822- try{
1823- int r = nonRedundantQIndeces.size() + i;
1824- int moI = redundantQIndeces[i].moI;
1825- int moJ = redundantQIndeces[i].moJ;
1826- if(moI == moJ){
1827- int rr = moI - (numberOcc-numberActiveOcc);
1828- q[r] = delta[rr];
1829- }
1830- else{
1831- q[r] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta)
1832- -this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
1833- }
1834- }
1835- catch(MolDSException ex){
1836-#pragma omp critical
1837- ompErrors << ex.what() << endl ;
1838- }
1839- }
1840- // Exception throwing for omp-region
1841- if(!ompErrors.str().empty()){
1842- throw MolDSException(ompErrors.str());
1843- }
1844- /*
1845- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1846- this->OutputLog(boost::format("q[%d] = %e\n") % i % q[i]);
1847- }
1848- for(int i=0; i<redundantQIndeces.size(); i++){
1849- int r = nonRedundantQIndeces.size() + i;
1850- this->OutputLog(boost::format("q[%d] = %e\n") % r % q[r]);
1851- }
1852- */
1853-}
1854-
1855-// see (40) and (45) in [PT_1996].
1856-// This method calculates "\Gamma_{NR} - K_{NR}" to solve (54) in [PT_1966]
1857-// Note taht K_{NR} is not calculated.
1858-void Mndo::CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR, const vector<MoIndexPair>& nonRedundantQIndeces) const{
1859- stringstream ompErrors;
1860-#pragma omp parallel for schedule(auto)
1861- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1862- try{
1863- int moI = nonRedundantQIndeces[i].moI;
1864- int moJ = nonRedundantQIndeces[i].moJ;
1865- for(int j=i; j<nonRedundantQIndeces.size(); j++){
1866- int moK = nonRedundantQIndeces[j].moI;
1867- int moL = nonRedundantQIndeces[j].moJ;
1868- gammaNRMinusKNR[i][j] = this->GetGammaNRElement(moI, moJ, moK, moL)
1869- -this->GetKNRElement(moI, moJ, moK, moL);
1870- }
1871- }
1872- catch(MolDSException ex){
1873-#pragma omp critical
1874- ompErrors << ex.what() << endl ;
1875- }
1876- }
1877- // Exception throwing for omp-region
1878- if(!ompErrors.str().empty()){
1879- throw MolDSException(ompErrors.str());
1880- }
1881-}
1882-
1883-// see (41), (42), and (46) in [PT_1996].
1884-// This method calculates "K_{R}^{\dagger} * Gamma_{R}" matrix, see (41), (42), and (46) to solve (54) in [PT_1996]
1885-// Note taht K_{R}^{\dager} is not calculated.
1886-void Mndo::CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
1887- const vector<MoIndexPair>& nonRedundantQIndeces,
1888- const vector<MoIndexPair>& redundantQIndeces) const{
1889- stringstream ompErrors;
1890-#pragma omp parallel for schedule(auto)
1891- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1892- try{
1893- int moI = nonRedundantQIndeces[i].moI;
1894- int moJ = nonRedundantQIndeces[i].moJ;
1895- for(int j=0; j<redundantQIndeces.size(); j++){
1896- int moK = redundantQIndeces[j].moI;
1897- int moL = redundantQIndeces[j].moJ;
1898- kRDagerGammaRInv[i][j] = this->GetKRDagerElement(moI, moJ, moK, moL)
1899- /this->GetGammaRElement(moK, moL, moK, moL);
1900- }
1901- }
1902- catch(MolDSException ex){
1903-#pragma omp critical
1904- ompErrors << ex.what() << endl ;
1905- }
1906- }
1907- // Exception throwing for omp-region
1908- if(!ompErrors.str().empty()){
1909- throw MolDSException(ompErrors.str());
1910- }
1911-}
1912-
1913-// right hand side of (54) in [PT_1996]
1914-void Mndo::CalcAuxiliaryVector(double* y,
1915- double const* q,
1916- double const* const* kRDagerGammaRInv,
1917- const vector<MoIndexPair>& nonRedundantQIndeces,
1918- const vector<MoIndexPair>& redundantQIndeces) const{
1919- MallocerFreer::GetInstance()->Initialize<double>(
1920- y,
1921- nonRedundantQIndeces.size());
1922- stringstream ompErrors;
1923-#pragma omp parallel for schedule(auto)
1924- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1925- try{
1926- int moI = nonRedundantQIndeces[i].moI;
1927- int moJ = nonRedundantQIndeces[i].moJ;
1928- y[i] += q[i]/this->GetNNRElement(moI, moJ, moI, moJ);
1929- for(int j=0; j<redundantQIndeces.size(); j++){
1930- int k = nonRedundantQIndeces.size() + j;
1931- y[i] += kRDagerGammaRInv[i][j]*q[k];
1932- }
1933- }
1934- catch(MolDSException ex){
1935-#pragma omp critical
1936- ompErrors << ex.what() << endl ;
1937- }
1938- }
1939- // Exception throwing for omp-region
1940- if(!ompErrors.str().empty()){
1941- throw MolDSException(ompErrors.str());
1942- }
1943-}
1944-
1945-void Mndo::TransposeFockMatrixMatrix(double** transposedFockMatrix) const{
1946- for(int i=0; i<this->molecule->GetTotalNumberAOs(); i++){
1947- for(int j=0; j<this->molecule->GetTotalNumberAOs(); j++){
1948- transposedFockMatrix[j][i] = this->fockMatrix[i][j];
1949- }
1950- }
1951-}
1952-
1953-// each element (mu, nu) of z matrix.
1954-// see (57) in [PT_1996]
1955-double Mndo::GetZMatrixForceElement(double const* y,
1956- double const* q,
1957- double const* const* transposedFockMatrix,
1958- const vector<MoIndexPair>& nonRedundantQIndeces,
1959- const vector<MoIndexPair>& redundantQIndeces,
1960- int mu,
1961- int nu) const{
1962- double value=0.0;
1963- for(int i=0; i<nonRedundantQIndeces.size(); i++){
1964- int moI = nonRedundantQIndeces[i].moI;
1965- int moJ = nonRedundantQIndeces[i].moJ;
1966- value += y[i]
1967- *transposedFockMatrix[mu][moI]
1968- *transposedFockMatrix[nu][moJ];
1969- }
1970- for(int i=0; i<redundantQIndeces.size(); i++){
1971- int j = nonRedundantQIndeces.size() + i;
1972- int moI = redundantQIndeces[i].moI;
1973- int moJ = redundantQIndeces[i].moJ;
1974- value += (q[j]/this->GetGammaRElement(moI, moJ, moI, moJ))
1975- *transposedFockMatrix[mu][moI]
1976- *transposedFockMatrix[nu][moJ];
1977- }
1978- return value;
1979-}
1980-
1981-void Mndo::CalcXiMatrices(double** xiOcc,
1982- double** xiVir,
1983- int exciteState,
1984- double const* const* transposedFockMatrix) const{
1985- int numberAOs = this->molecule->GetTotalNumberAOs();
1986- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
1987- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
1988- int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
1989- MallocerFreer::GetInstance()->Initialize<double>(
1990- xiOcc, numberActiveOcc, numberAOs);
1991- MallocerFreer::GetInstance()->Initialize<double>(
1992- xiVir, numberActiveVir, numberAOs);
1993- stringstream ompErrors;
1994- // xiOcc
1995-#pragma omp parallel for schedule(auto)
1996- for(int p=0; p<numberActiveOcc; p++){
1997- try{
1998- for(int mu=0; mu<numberAOs; mu++){
1999- for(int a=0; a<numberActiveVir; a++){
2000- int moA = numberOcc + a;
2001- int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(p,a);
2002- xiOcc[p][mu] += this->matrixCIS[exciteState][slaterDeterminantIndex]
2003- *transposedFockMatrix[mu][moA];
2004- }
2005- }
2006- }
2007- catch(MolDSException ex){
2008-#pragma omp critical
2009- ompErrors << ex.what() << endl ;
2010- }
2011- }
2012- // Exception throwing for omp-region
2013- if(!ompErrors.str().empty()){
2014- throw MolDSException(ompErrors.str());
2015- }
2016- // xiVir
2017-#pragma omp parallel for schedule(auto)
2018- for(int p=0; p<numberActiveVir; p++){
2019- try{
2020- for(int mu=0; mu<numberAOs; mu++){
2021- for(int i=0; i<numberActiveOcc; i++){
2022- int moI = numberOcc - (i+1);
2023- int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,p);
2024- xiVir[p][mu] += this->matrixCIS[exciteState][slaterDeterminantIndex]
2025- *transposedFockMatrix[mu][moI];
2026- }
2027- }
2028- }
2029- catch(MolDSException ex){
2030-#pragma omp critical
2031- ompErrors << ex.what() << endl ;
2032- }
2033- }
2034- // Exception throwing for omp-region
2035- if(!ompErrors.str().empty()){
2036- throw MolDSException(ompErrors.str());
2037- }
2038-}
2039-
2040941 void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
2041942 double****** diatomicOverlapAOs2ndDerivs,
2042943 double******* diatomicTwoElecTwoCore1stDerivs,
@@ -3148,182 +2049,6 @@ void Mndo::FreeTempMatricesSolveCPHF(double*** matrixCPHF,
31482049 MallocerFreer::GetInstance()->Free<double>(matrixCPHF, dimensionCPHF, dimensionCPHF);
31492050 }
31502051
3151-// see [PT_1996, PT_1997]
3152-void Mndo::CalcZMatrixForce(const vector<int>& elecStates){
3153-#ifdef MOLDS_DBG
3154- if(this->etaMatrixForce == NULL){
3155- throw MolDSException(this->errorMessageCalcZMatrixForceEtaNull);
3156- }
3157-#endif
3158- this->CheckZMatrixForce(elecStates);
3159-
3160- // creat MO-index-pair for Q variables.
3161- vector<MoIndexPair> nonRedundantQIndeces;
3162- vector<MoIndexPair> redundantQIndeces;
3163- this->CalcActiveSetVariablesQ(&nonRedundantQIndeces,
3164- &redundantQIndeces,
3165- Parameters::GetInstance()->GetActiveOccCIS(),
3166- Parameters::GetInstance()->GetActiveVirCIS());
3167-
3168- // malloc temporary arrays
3169- double* delta = NULL; // Delta matrix, see (9) in [PT_1997]
3170- double* q = NULL; //// Q-vector in (19) in [PT_1997]
3171- double** gammaNRMinusKNR = NULL; // Gmamma_{NR} - K_{NR} matrix, see (40) and (45) to slove (54) in [PT_1996]
3172- double** kRDagerGammaRInv = NULL; // K_{R}^{\dagger} * Gamma_{R} matrix, see (41), (42), and (46) to solve (54) in [PT_1996]
3173- double* y = NULL; // y-vector in (54) in [PT_1996]
3174- double** transposedFockMatrix = NULL; // transposed Fock matrix
3175- double** xiOcc = NULL;
3176- double** xiVir = NULL;
3177- try{
3178- this->MallocTempMatrixForZMatrix(&delta,
3179- &q,
3180- &gammaNRMinusKNR,
3181- &kRDagerGammaRInv,
3182- &y,
3183- &transposedFockMatrix,
3184- &xiOcc,
3185- &xiVir,
3186- nonRedundantQIndeces.size(),
3187- redundantQIndeces.size());
3188- this->TransposeFockMatrixMatrix(transposedFockMatrix);
3189- this->CalcGammaNRMinusKNRMatrix(gammaNRMinusKNR, nonRedundantQIndeces);
3190- this->CalcKRDagerGammaRInvMatrix(kRDagerGammaRInv, nonRedundantQIndeces,redundantQIndeces);
3191- int groundState=0;
3192- for(int n=0; n<elecStates.size(); n++){
3193- if(groundState < elecStates[n]){
3194- int exciteState = elecStates[n]-1;
3195- this->CalcDeltaVector(delta, exciteState);
3196- this->CalcXiMatrices(xiOcc, xiVir, exciteState, transposedFockMatrix);
3197- this->CalcQVector(q,
3198- delta,
3199- xiOcc,
3200- xiVir,
3201- this->etaMatrixForce[n],
3202- nonRedundantQIndeces,
3203- redundantQIndeces);
3204- this->CalcAuxiliaryVector(y, q, kRDagerGammaRInv, nonRedundantQIndeces, redundantQIndeces);
3205- // solve (54) in [PT_1996]
3206- MolDS_wrappers::Lapack::GetInstance()->Dsysv(gammaNRMinusKNR,
3207- y,
3208- nonRedundantQIndeces.size());
3209- // calculate each element of Z matrix.
3210- stringstream ompErrors;
3211-#pragma omp parallel for schedule(auto)
3212- for(int mu=0; mu<this->molecule->GetTotalNumberAOs(); mu++){
3213- try{
3214- for(int nu=0; nu<this->molecule->GetTotalNumberAOs(); nu++){
3215- this->zMatrixForce[n][mu][nu] = this->GetZMatrixForceElement(
3216- y,
3217- q,
3218- transposedFockMatrix,
3219- nonRedundantQIndeces,
3220- redundantQIndeces,
3221- mu,
3222- nu);
3223- }
3224- }
3225- catch(MolDSException ex){
3226-#pragma omp critical
3227- ompErrors << ex.what() << endl ;
3228- }
3229- }
3230- // Exception throwing for omp-region
3231- if(!ompErrors.str().empty()){
3232- throw MolDSException(ompErrors.str());
3233- }
3234-
3235- }
3236- }
3237- }
3238- catch(MolDSException ex){
3239- this->FreeTempMatrixForZMatrix(&delta,
3240- &q,
3241- &gammaNRMinusKNR,
3242- &kRDagerGammaRInv,
3243- &y,
3244- &transposedFockMatrix,
3245- &xiOcc,
3246- &xiVir,
3247- nonRedundantQIndeces.size(),
3248- redundantQIndeces.size());
3249- throw ex;
3250- }
3251- this->FreeTempMatrixForZMatrix(&delta,
3252- &q,
3253- &gammaNRMinusKNR,
3254- &kRDagerGammaRInv,
3255- &y,
3256- &transposedFockMatrix,
3257- &xiOcc,
3258- &xiVir,
3259- nonRedundantQIndeces.size(),
3260- redundantQIndeces.size());
3261-}
3262-
3263-void Mndo::CalcEtaMatrixForce(const vector<int>& elecStates){
3264- this->CheckEtaMatrixForce(elecStates);
3265- int numberAOs = this->molecule->GetTotalNumberAOs();
3266- int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3267- int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
3268- int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
3269- int groundState = 0;
3270- double** transposedFockMatrix = NULL; // transposed Fock matrix
3271- try{
3272- MallocerFreer::GetInstance()->Malloc<double>(&transposedFockMatrix,
3273- numberAOs,
3274- numberAOs);
3275- this->TransposeFockMatrixMatrix(transposedFockMatrix);
3276- for(int n=0; n<elecStates.size(); n++){
3277- if(groundState < elecStates[n]){
3278- int exciteState = elecStates[n]-1;
3279-
3280- // calc each element
3281- stringstream ompErrors;
3282-#pragma omp parallel for schedule(auto)
3283- for(int mu=0; mu<numberAOs; mu++){
3284- try{
3285- for(int nu=0; nu<numberAOs; nu++){
3286- for(int i=0; i<numberActiveOcc; i++){
3287- int moI = numberOcc-(i+1);
3288- for(int a=0; a<numberActiveVir; a++){
3289- int moA = numberOcc+a;
3290- int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,a);
3291- this->etaMatrixForce[n][mu][nu]
3292- += this->matrixCIS[exciteState][slaterDeterminantIndex]
3293- *transposedFockMatrix[mu][moI]
3294- *transposedFockMatrix[nu][moA];
3295- }
3296- }
3297- }
3298- }
3299- catch(MolDSException ex){
3300-#pragma omp critical
3301- ompErrors << ex.what() << endl ;
3302- }
3303- }
3304- // Exception throwing for omp-region
3305- if(!ompErrors.str().empty()){
3306- throw MolDSException(ompErrors.str());
3307- }
3308-
3309- }
3310- }
3311- }
3312- catch(MolDSException ex){
3313- MallocerFreer::GetInstance()->Free<double>(&transposedFockMatrix,numberAOs,numberAOs);
3314- throw ex;
3315- }
3316- MallocerFreer::GetInstance()->Free<double>(&transposedFockMatrix,numberAOs, numberAOs);
3317-}
3318-
3319-bool Mndo::RequiresExcitedStatesForce(const vector<int>& elecStates) const{
3320- bool requires = true;
3321- if(elecStates.size()==1 && elecStates[0]==0){
3322- requires = false;
3323- }
3324- return requires;
3325-}
3326-
33272052 void Mndo::CalcForceSCFElecCoreAttractionPart(double* force,
33282053 int indexAtomA,
33292054 int indexAtomB,
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -43,7 +43,6 @@ protected:
4343 std::string errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms;
4444 std::string errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix;
4545 std::string errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix;
46- std::string errorMessageCalcZMatrixForceEtaNull;
4746 virtual void SetMessages();
4847 virtual void SetEnableAtomTypes();
4948 virtual void CalcSCFProperties();
@@ -107,12 +106,7 @@ private:
107106 std::string errorMessageMultipoleB;
108107 std::string messageHeatsFormation;
109108 std::string messageHeatsFormationTitle;
110- struct MoIndexPair{int moI; int moJ; bool isMoICIMO; bool isMoJCIMO;};
111109 double heatsFormation;
112- int zMatrixForceElecStatesNum;
113- int etaMatrixForceElecStatesNum;
114- double*** zMatrixForce;
115- double*** etaMatrixForce;
116110 double GetAuxiliaryDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
117111 const MolDS_base_atoms::Atom& atomB,
118112 double distanceAB) const;
@@ -125,80 +119,8 @@ private:
125119 double distanceAB,
126120 MolDS_base::CartesianType axisA1,
127121 MolDS_base::CartesianType axisA2) const;
128- double GetGammaNRElement(int moI, int moJ, int moK, int moL) const;
129- double GetGammaRElement(int moI, int moJ, int moK, int moL) const;
130- double GetNNRElement(int moI, int moJ, int moK, int moL) const;
131- double GetNRElement(int moI, int moJ, int moK, int moL) const;
132- double GetKNRElement(int moI, int moJ, int moK, int moL) const;
133- double GetKRElement(int moI, int moJ, int moK, int moL) const;
134- double GetKRDagerElement(int moI, int moJ, int moK, int moL) const;
135- double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
136- void MallocTempMatrixForZMatrix(double** delta,
137- double** q,
138- double*** gammaNRMinusKNR,
139- double*** kRDag,
140- double** y,
141- double*** transposedFockMatrix,
142- double*** xiOcc,
143- double*** xiVir,
144- int sizeQNR,
145- int sizeQR) const;
146- void FreeTempMatrixForZMatrix(double** delta,
147- double** q,
148- double*** gammaNRMinusKNR,
149- double*** kRDag,
150- double** y,
151- double*** transposedFockMatrix,
152- double*** xiOcc,
153- double*** xiVir,
154- int sizeQNR,
155- int sizeQR) const;
156- void CalcDeltaVector(double* delta, int exciteState) const;
157- double GetSmallQElement(int moI,
158- int moP,
159- double const* const* xiOcc,
160- double const* const* xiVir,
161- double const* const* eta) const;
162- void CalcQVector(double* q,
163- double const* delta,
164- double const* const* xiOcc,
165- double const* const* xiVir,
166- double const* const* eta,
167- const std::vector<MoIndexPair>& nonRedundantQIndeces,
168- const std::vector<MoIndexPair>& redundantQIndeces) const;
169- void TransposeFockMatrixMatrix(double** transposedFockMatrix) const;
170- void CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR,
171- const std::vector<MoIndexPair>& nonRedundantQIndeces) const;
172- void CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
173- const std::vector<MoIndexPair>& nonRedundantQIndeces,
174- const std::vector<MoIndexPair>& redundantQIndeces) const;
175- void CalcAuxiliaryVector(double* y,
176- double const* q,
177- double const* const* kRDagerGammaRInv,
178- const std::vector<MoIndexPair>& nonRedundantQIndeces,
179- const std::vector<MoIndexPair>& redundantQIndeces) const;
180- void CalcXiMatrices(double** xiOcc,
181- double** xiVir,
182- int exciteState,
183- double const* const* transposedFockMatrix) const;
184- double GetZMatrixForceElement(double const* y,
185- double const* q,
186- double const* const* transposedFockMatrix,
187- const std::vector<MoIndexPair>& nonRedundantQIndeces,
188- const std::vector<MoIndexPair>& redundantQIndeces,
189- int mu,
190- int nu) const;
191- void CheckZMatrixForce(const std::vector<int>& elecStates);
192- void CheckEtaMatrixForce(const std::vector<int>& elecStates);
193- void CalcZMatrixForce(const std::vector<int>& elecStates);
194- void CalcEtaMatrixForce(const std::vector<int>& elecStates);
195- bool RequiresExcitedStatesForce(const std::vector<int>& elecStates) const;
196122 double GetCISCoefficientMOEnergy(int k, int l, int r, int numberActiveVir) const;
197123 double GetCISCoefficientTwoElecIntegral(int k, int l, int p, int q, int r, int s, int numberActiveVir) const;
198- void CalcActiveSetVariablesQ(std::vector<MoIndexPair>* nonRedundantQIndeces,
199- std::vector<MoIndexPair>* redundantQIndeces,
200- int numberActiveOcc,
201- int numberActiveVir) const;
202124 void CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const;
203125 double GetHessianElementSameAtomsSCF(int indexAtomA,
204126 MolDS_base::CartesianType axisA1,
--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -67,6 +67,11 @@ ZindoS::ZindoS() : MolDS_cndo::Cndo2(){
6767 this->SetMessages();
6868 this->SetEnableAtomTypes();
6969
70+ this->zMatrixForceElecStatesNum = 0;
71+ this->etaMatrixForceElecStatesNum = 0;
72+ this->zMatrixForce = NULL;
73+ this->etaMatrixForce = NULL;
74+
7075 //private variables
7176 this->matrixForceElecStatesNum = 0;
7277 this->nishimotoMatagaParamA = 1.2;
@@ -88,6 +93,14 @@ ZindoS::~ZindoS(){
8893 this->matrixForceElecStatesNum,
8994 this->molecule->GetNumberAtoms(),
9095 CartesianType_end);
96+ MallocerFreer::GetInstance()->Free<double>(&this->zMatrixForce,
97+ this->zMatrixForceElecStatesNum,
98+ this->molecule->GetTotalNumberAOs(),
99+ this->molecule->GetTotalNumberAOs());
100+ MallocerFreer::GetInstance()->Free<double>(&this->etaMatrixForce,
101+ this->etaMatrixForceElecStatesNum,
102+ this->molecule->GetTotalNumberAOs(),
103+ this->molecule->GetTotalNumberAOs());
91104 if(Parameters::GetInstance()->RequiresMullikenCIS()){
92105 vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
93106 MallocerFreer::GetInstance()->Free<double>(&this->orbitalElectronPopulationCIS,
@@ -138,6 +151,8 @@ void ZindoS::SetMessages(){
138151 = "Error in zindo::ZindoS::CalcElectronicTransitionDipoleMoment: Bad eigen state is set to calculate the transition dipole moment. Note taht state=0 means the ground state and other state = i means the i-th excited state in below.\n";
139152 this->errorMessageCalcFrequenciesNormalModesBadTheory
140153 = "Error in zindo::ZindoS::CalcFrequenciesNormalModesBadTheory: ZINDO/S is not supported for frequency (normal mode) analysis.\n";
154+ this->errorMessageCalcZMatrixForceEtaNull
155+ = "Error in zindo::ZindoS::CalcZMatrixForce: Nndo::etaMatrixForce is NULL. Call Mndo::CalcEtaMatrixForce before calling Mndo::CalcZMatrixForce.\n";
141156 this->messageSCFMetConvergence = "\n\n\n\t\tZINDO/S-SCF met convergence criterion(^^b\n\n\n";
142157 this->messageStartSCF = "********** START: ZINDO/S-SCF **********\n";
143158 this->messageDoneSCF = "********** DONE: ZINDO/S-SCF **********\n\n\n";
@@ -2547,6 +2562,1261 @@ int ZindoS::GetActiveVirIndex(const MolDS_base::Molecule& molecule, int matrixCI
25472562 +(matrixCISIndex%Parameters::GetInstance()->GetActiveVirCIS());
25482563 }
25492564
2565+bool ZindoS::RequiresExcitedStatesForce(const vector<int>& elecStates) const{
2566+ bool requires = true;
2567+ if(elecStates.size()==1 && elecStates[0]==0){
2568+ requires = false;
2569+ }
2570+ return requires;
2571+}
2572+
2573+void ZindoS::CheckZMatrixForce(const vector<int>& elecStates){
2574+ // malloc or initialize Z matrix
2575+ if(this->zMatrixForce == NULL){
2576+ MallocerFreer::GetInstance()->Malloc<double>(&this->zMatrixForce,
2577+ elecStates.size(),
2578+ this->molecule->GetTotalNumberAOs(),
2579+ this->molecule->GetTotalNumberAOs());
2580+ this->zMatrixForceElecStatesNum = elecStates.size();
2581+ }
2582+ else{
2583+ MallocerFreer::GetInstance()->
2584+ Initialize<double>(this->zMatrixForce,
2585+ elecStates.size(),
2586+ this->molecule->GetTotalNumberAOs(),
2587+ this->molecule->GetTotalNumberAOs());
2588+ }
2589+}
2590+
2591+void ZindoS::CheckEtaMatrixForce(const vector<int>& elecStates){
2592+ // malloc or initialize eta matrix
2593+ if(this->etaMatrixForce == NULL){
2594+ MallocerFreer::GetInstance()->Malloc<double>(&this->etaMatrixForce,
2595+ elecStates.size(),
2596+ this->molecule->GetTotalNumberAOs(),
2597+ this->molecule->GetTotalNumberAOs());
2598+ this->etaMatrixForceElecStatesNum = elecStates.size();
2599+ }
2600+ else{
2601+ MallocerFreer::GetInstance()->
2602+ Initialize<double>(this->etaMatrixForce,
2603+ elecStates.size(),
2604+ this->molecule->GetTotalNumberAOs(),
2605+ this->molecule->GetTotalNumberAOs());
2606+ }
2607+}
2608+
2609+void ZindoS::CalcEtaMatrixForce(const vector<int>& elecStates){
2610+ this->CheckEtaMatrixForce(elecStates);
2611+ int numberAOs = this->molecule->GetTotalNumberAOs();
2612+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
2613+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
2614+ int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
2615+ int groundState = 0;
2616+ double** transposedFockMatrix = NULL; // transposed Fock matrix
2617+ try{
2618+ MallocerFreer::GetInstance()->Malloc<double>(&transposedFockMatrix,
2619+ numberAOs,
2620+ numberAOs);
2621+ this->TransposeFockMatrixMatrix(transposedFockMatrix);
2622+ for(int n=0; n<elecStates.size(); n++){
2623+ if(groundState < elecStates[n]){
2624+ int exciteState = elecStates[n]-1;
2625+
2626+ // calc each element
2627+ stringstream ompErrors;
2628+#pragma omp parallel for schedule(auto)
2629+ for(int mu=0; mu<numberAOs; mu++){
2630+ try{
2631+ for(int nu=0; nu<numberAOs; nu++){
2632+ for(int i=0; i<numberActiveOcc; i++){
2633+ int moI = numberOcc-(i+1);
2634+ for(int a=0; a<numberActiveVir; a++){
2635+ int moA = numberOcc+a;
2636+ int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,a);
2637+ this->etaMatrixForce[n][mu][nu]
2638+ += this->matrixCIS[exciteState][slaterDeterminantIndex]
2639+ *transposedFockMatrix[mu][moI]
2640+ *transposedFockMatrix[nu][moA];
2641+ }
2642+ }
2643+ }
2644+ }
2645+ catch(MolDSException ex){
2646+#pragma omp critical
2647+ ompErrors << ex.what() << endl ;
2648+ }
2649+ }
2650+ // Exception throwing for omp-region
2651+ if(!ompErrors.str().empty()){
2652+ throw MolDSException(ompErrors.str());
2653+ }
2654+
2655+ }
2656+ }
2657+ }
2658+ catch(MolDSException ex){
2659+ MallocerFreer::GetInstance()->Free<double>(&transposedFockMatrix,numberAOs,numberAOs);
2660+ throw ex;
2661+ }
2662+ MallocerFreer::GetInstance()->Free<double>(&transposedFockMatrix,numberAOs, numberAOs);
2663+}
2664+
2665+// see [PT_1996, PT_1997]
2666+void ZindoS::CalcZMatrixForce(const vector<int>& elecStates){
2667+#ifdef MOLDS_DBG
2668+ if(this->etaMatrixForce == NULL){
2669+ throw MolDSException(this->errorMessageCalcZMatrixForceEtaNull);
2670+ }
2671+#endif
2672+ this->CheckZMatrixForce(elecStates);
2673+
2674+ // creat MO-index-pair for Q variables.
2675+ vector<MoIndexPair> nonRedundantQIndeces;
2676+ vector<MoIndexPair> redundantQIndeces;
2677+ this->CalcActiveSetVariablesQ(&nonRedundantQIndeces,
2678+ &redundantQIndeces,
2679+ Parameters::GetInstance()->GetActiveOccCIS(),
2680+ Parameters::GetInstance()->GetActiveVirCIS());
2681+
2682+ // malloc temporary arrays
2683+ double* delta = NULL; // Delta matrix, see (9) in [PT_1997]
2684+ double* q = NULL; //// Q-vector in (19) in [PT_1997]
2685+ double** gammaNRMinusKNR = NULL; // Gmamma_{NR} - K_{NR} matrix, see (40) and (45) to slove (54) in [PT_1996]
2686+ double** kRDagerGammaRInv = NULL; // K_{R}^{\dagger} * Gamma_{R} matrix, see (41), (42), and (46) to solve (54) in [PT_1996]
2687+ double* y = NULL; // y-vector in (54) in [PT_1996]
2688+ double** transposedFockMatrix = NULL; // transposed Fock matrix
2689+ double** xiOcc = NULL;
2690+ double** xiVir = NULL;
2691+ try{
2692+ this->MallocTempMatrixForZMatrix(&delta,
2693+ &q,
2694+ &gammaNRMinusKNR,
2695+ &kRDagerGammaRInv,
2696+ &y,
2697+ &transposedFockMatrix,
2698+ &xiOcc,
2699+ &xiVir,
2700+ nonRedundantQIndeces.size(),
2701+ redundantQIndeces.size());
2702+ this->TransposeFockMatrixMatrix(transposedFockMatrix);
2703+ this->CalcGammaNRMinusKNRMatrix(gammaNRMinusKNR, nonRedundantQIndeces);
2704+ this->CalcKRDagerGammaRInvMatrix(kRDagerGammaRInv, nonRedundantQIndeces,redundantQIndeces);
2705+ int groundState=0;
2706+ 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;
2725+#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 ;
2742+ }
2743+ }
2744+ // Exception throwing for omp-region
2745+ if(!ompErrors.str().empty()){
2746+ throw MolDSException(ompErrors.str());
2747+ }
2748+
2749+ }
2750+ }
2751+ }
2752+ catch(MolDSException ex){
2753+ this->FreeTempMatrixForZMatrix(&delta,
2754+ &q,
2755+ &gammaNRMinusKNR,
2756+ &kRDagerGammaRInv,
2757+ &y,
2758+ &transposedFockMatrix,
2759+ &xiOcc,
2760+ &xiVir,
2761+ nonRedundantQIndeces.size(),
2762+ redundantQIndeces.size());
2763+ throw ex;
2764+ }
2765+ this->FreeTempMatrixForZMatrix(&delta,
2766+ &q,
2767+ &gammaNRMinusKNR,
2768+ &kRDagerGammaRInv,
2769+ &y,
2770+ &transposedFockMatrix,
2771+ &xiOcc,
2772+ &xiVir,
2773+ nonRedundantQIndeces.size(),
2774+ redundantQIndeces.size());
2775+}
2776+
2777+// each element (mu, nu) of z matrix.
2778+// see (57) in [PT_1996]
2779+double ZindoS::GetZMatrixForceElement(double const* y,
2780+ double const* q,
2781+ double const* const* transposedFockMatrix,
2782+ const vector<MoIndexPair>& nonRedundantQIndeces,
2783+ const vector<MoIndexPair>& redundantQIndeces,
2784+ int mu,
2785+ int nu) const{
2786+ double value=0.0;
2787+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
2788+ int moI = nonRedundantQIndeces[i].moI;
2789+ int moJ = nonRedundantQIndeces[i].moJ;
2790+ value += y[i]
2791+ *transposedFockMatrix[mu][moI]
2792+ *transposedFockMatrix[nu][moJ];
2793+ }
2794+ for(int i=0; i<redundantQIndeces.size(); i++){
2795+ int j = nonRedundantQIndeces.size() + i;
2796+ int moI = redundantQIndeces[i].moI;
2797+ int moJ = redundantQIndeces[i].moJ;
2798+ value += (q[j]/this->GetGammaRElement(moI, moJ, moI, moJ))
2799+ *transposedFockMatrix[mu][moI]
2800+ *transposedFockMatrix[nu][moJ];
2801+ }
2802+ return value;
2803+}
2804+
2805+void ZindoS::MallocTempMatrixForZMatrix(double** delta,
2806+ double** q,
2807+ double*** gammaNRMinusKNR,
2808+ double*** kRDag,
2809+ double** y,
2810+ double*** transposedFockMatrix,
2811+ double*** xiOcc,
2812+ double*** xiVir,
2813+ int sizeQNR,
2814+ int sizeQR) const{
2815+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
2816+ int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
2817+ int numberActiveMO = numberActiveOcc + numberActiveVir;
2818+ int numberAOs = this->molecule->GetTotalNumberAOs();
2819+ MallocerFreer::GetInstance()->Malloc<double>(delta, numberActiveMO);
2820+ MallocerFreer::GetInstance()->Malloc<double>(q, sizeQNR+sizeQR);
2821+ MallocerFreer::GetInstance()->Malloc<double>(gammaNRMinusKNR, sizeQNR, sizeQNR);
2822+ MallocerFreer::GetInstance()->Malloc<double>(kRDag, sizeQNR, sizeQR);
2823+ MallocerFreer::GetInstance()->Malloc<double>(y, sizeQNR);
2824+ MallocerFreer::GetInstance()->Malloc<double>(transposedFockMatrix,
2825+ numberAOs,
2826+ numberAOs);
2827+ MallocerFreer::GetInstance()->Malloc<double>(xiOcc, numberActiveOcc,numberAOs);
2828+ MallocerFreer::GetInstance()->Malloc<double>(xiVir,numberActiveVir,numberAOs);
2829+}
2830+
2831+void ZindoS::FreeTempMatrixForZMatrix(double** delta,
2832+ double** q,
2833+ double*** gammaNRMinusKNR,
2834+ double*** kRDag,
2835+ double** y,
2836+ double*** transposedFockMatrix,
2837+ double*** xiOcc,
2838+ double*** xiVir,
2839+ int sizeQNR,
2840+ int sizeQR) const{
2841+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
2842+ int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
2843+ int numberActiveMO = numberActiveOcc + numberActiveVir;
2844+ int numberAOs = this->molecule->GetTotalNumberAOs();
2845+ MallocerFreer::GetInstance()->Free<double>(delta, numberActiveMO);
2846+ MallocerFreer::GetInstance()->Free<double>(q, sizeQNR+sizeQR);
2847+ MallocerFreer::GetInstance()->Free<double>(gammaNRMinusKNR, sizeQNR, sizeQNR);
2848+ MallocerFreer::GetInstance()->Free<double>(kRDag, sizeQNR, sizeQR);
2849+ MallocerFreer::GetInstance()->Free<double>(y, sizeQNR);
2850+ MallocerFreer::GetInstance()->Free<double>(transposedFockMatrix, numberAOs, numberAOs);
2851+ MallocerFreer::GetInstance()->Free<double>(xiOcc, numberActiveOcc, numberAOs);
2852+ MallocerFreer::GetInstance()->Free<double>(xiVir, numberActiveVir, numberAOs);
2853+}
2854+
2855+// see (9) in [PT_1997]
2856+void ZindoS::CalcDeltaVector(double* delta, int exciteState) const{
2857+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
2858+ int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
2859+ int numberActiveMO = numberActiveOcc + numberActiveVir;
2860+ MallocerFreer::GetInstance()->Initialize<double>(delta, numberActiveMO);
2861+ stringstream ompErrors;
2862+#pragma omp parallel for schedule(auto)
2863+ for(int r=0; r<numberActiveMO; r++){
2864+ try{
2865+ double value = 0.0;
2866+ if(r<numberActiveOcc){
2867+ // r is active occupied MO
2868+ int rr=numberActiveOcc-(r+1);
2869+ for(int a=0; a<numberActiveVir; a++){
2870+ int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(rr,a);
2871+ value -= pow(this->matrixCIS[exciteState][slaterDeterminantIndex],2.0);
2872+ }
2873+ }
2874+ else{
2875+ // r is active virtual MO
2876+ int rr=r-numberActiveOcc;
2877+ for(int i=0; i<numberActiveOcc; i++){
2878+ int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,rr);
2879+ value += pow(this->matrixCIS[exciteState][slaterDeterminantIndex],2.0);
2880+ }
2881+ }
2882+ delta[r] = value;
2883+ }
2884+ catch(MolDSException ex){
2885+#pragma omp critical
2886+ ompErrors << ex.what() << endl ;
2887+ }
2888+ }
2889+ // Exception throwing for omp-region
2890+ if(!ompErrors.str().empty()){
2891+ throw MolDSException(ompErrors.str());
2892+ }
2893+}
2894+
2895+// see variable Q-vector in [PT_1996, PT_1997]
2896+void ZindoS::CalcActiveSetVariablesQ(vector<MoIndexPair>* nonRedundantQIndeces,
2897+ vector<MoIndexPair>* redundantQIndeces,
2898+ int numberActiveOcc,
2899+ int numberActiveVir) const{
2900+ int numberAOs = this->molecule->GetTotalNumberAOs();
2901+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
2902+ for(int moI=0; moI<numberOcc; moI++){
2903+ bool isMoICIMO = numberOcc-numberActiveOcc<=moI ? true : false;
2904+ for(int moJ=numberOcc; moJ<numberAOs; moJ++){
2905+ bool isMoJCIMO = moJ<numberOcc+numberActiveVir ? true : false;
2906+ MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
2907+ nonRedundantQIndeces->push_back(moIndexPair);
2908+ }
2909+ }
2910+ for(int moI=numberOcc-numberActiveOcc; moI<numberOcc; moI++){
2911+ bool isMoICIMO = true;
2912+ for(int moJ=moI; moJ<numberOcc; moJ++){
2913+ bool isMoJCIMO = true;
2914+ MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
2915+ redundantQIndeces->push_back(moIndexPair);
2916+ }
2917+ }
2918+ for(int moI=numberOcc; moI<numberOcc+numberActiveVir; moI++){
2919+ bool isMoICIMO = true;
2920+ for(int moJ=moI; moJ<numberOcc+numberActiveVir; moJ++){
2921+ bool isMoJCIMO = true;
2922+ MoIndexPair moIndexPair = {moI, moJ, isMoICIMO, isMoJCIMO};
2923+ redundantQIndeces->push_back(moIndexPair);
2924+ }
2925+ }
2926+}
2927+
2928+// see (20) - (23) in [PT_1997]
2929+void ZindoS::CalcQVector(double* q,
2930+ double const* delta,
2931+ double const* const* xiOcc,
2932+ double const* const* xiVir,
2933+ double const* const* eta,
2934+ const vector<MoIndexPair>& nonRedundantQIndeces,
2935+ const vector<MoIndexPair>& redundantQIndeces) const{
2936+ MallocerFreer::GetInstance()->Initialize<double>(
2937+ q,
2938+ nonRedundantQIndeces.size()+redundantQIndeces.size());
2939+
2940+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
2941+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
2942+ stringstream ompErrors;
2943+#pragma omp parallel for schedule(auto)
2944+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
2945+ try{
2946+ int moI = nonRedundantQIndeces[i].moI;
2947+ int moJ = nonRedundantQIndeces[i].moJ;
2948+ bool isMoICIMO = nonRedundantQIndeces[i].isMoICIMO;
2949+ bool isMoJCIMO = nonRedundantQIndeces[i].isMoJCIMO;
2950+ if(!isMoICIMO && isMoJCIMO){
2951+ q[i] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta);
2952+ }
2953+ else if(isMoICIMO && !isMoJCIMO){
2954+ q[i] = -1.0*this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
2955+ }
2956+ else if(isMoICIMO && isMoJCIMO){
2957+ q[i] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta)
2958+ -this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
2959+ }
2960+ else{
2961+ q[i] = 0.0;
2962+ }
2963+ }
2964+ catch(MolDSException ex){
2965+#pragma omp critical
2966+ ompErrors << ex.what() << endl ;
2967+ }
2968+ }
2969+ // Exception throwing for omp-region
2970+ if(!ompErrors.str().empty()){
2971+ throw MolDSException(ompErrors.str());
2972+ }
2973+#pragma omp parallel for schedule(auto)
2974+ for(int i=0; i<redundantQIndeces.size(); i++){
2975+ try{
2976+ int r = nonRedundantQIndeces.size() + i;
2977+ int moI = redundantQIndeces[i].moI;
2978+ int moJ = redundantQIndeces[i].moJ;
2979+ if(moI == moJ){
2980+ int rr = moI - (numberOcc-numberActiveOcc);
2981+ q[r] = delta[rr];
2982+ }
2983+ else{
2984+ q[r] = this->GetSmallQElement(moI, moJ, xiOcc, xiVir, eta)
2985+ -this->GetSmallQElement(moJ, moI, xiOcc, xiVir, eta);
2986+ }
2987+ }
2988+ catch(MolDSException ex){
2989+#pragma omp critical
2990+ ompErrors << ex.what() << endl ;
2991+ }
2992+ }
2993+ // Exception throwing for omp-region
2994+ if(!ompErrors.str().empty()){
2995+ throw MolDSException(ompErrors.str());
2996+ }
2997+ /*
2998+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
2999+ this->OutputLog(boost::format("q[%d] = %e\n") % i % q[i]);
3000+ }
3001+ for(int i=0; i<redundantQIndeces.size(); i++){
3002+ int r = nonRedundantQIndeces.size() + i;
3003+ this->OutputLog(boost::format("q[%d] = %e\n") % r % q[r]);
3004+ }
3005+ */
3006+}
3007+
3008+// see (18) in [PT_1977]
3009+double ZindoS::GetSmallQElement(int moI,
3010+ int moP,
3011+ double const* const* xiOcc,
3012+ double const* const* xiVir,
3013+ double const* const* eta) const{
3014+ double value = 0.0;
3015+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3016+ bool isMoPOcc = moP<numberOcc ? true : false;
3017+
3018+ for(int A=0; A<molecule->GetNumberAtoms(); A++){
3019+ const Atom& atomA = *molecule->GetAtom(A);
3020+ int firstAOIndexA = atomA.GetFirstAOIndex();
3021+ int lastAOIndexA = atomA.GetLastAOIndex();
3022+
3023+ for(int B=A; B<molecule->GetNumberAtoms(); B++){
3024+ const Atom& atomB = *molecule->GetAtom(B);
3025+ int firstAOIndexB = atomB.GetFirstAOIndex();
3026+ int lastAOIndexB = atomB.GetLastAOIndex();
3027+
3028+ if(A!=B){
3029+ 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+ }
3140+ }
3141+ }
3142+ }
3143+ }
3144+ else{
3145+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3146+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3147+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3148+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3149+ double twoElecInt = 0.0;
3150+ if(mu==nu && lambda==sigma){
3151+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3152+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3153+ twoElecInt = this->GetCoulombInt(orbitalMu,
3154+ orbitalLambda,
3155+ atomA);
3156+ }
3157+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3158+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3159+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
3160+ twoElecInt = this->GetExchangeInt(orbitalMu,
3161+ orbitalNu,
3162+ atomA);
3163+ }
3164+ else{
3165+ twoElecInt = 0.0;
3166+ }
3167+
3168+ double temp = 0.0;
3169+ if(isMoPOcc){
3170+ int p = numberOcc - (moP+1);
3171+ temp = 4.0*xiOcc[p][nu]*eta[lambda][sigma]
3172+ -1.0*xiOcc[p][lambda]*eta[nu][sigma]
3173+ -1.0*xiOcc[p][sigma]*eta[nu][lambda];
3174+ }
3175+ else{
3176+ int p = moP - numberOcc;
3177+ temp = 4.0*xiVir[p][nu]*eta[lambda][sigma]
3178+ -1.0*xiVir[p][lambda]*eta[sigma][nu]
3179+ -1.0*xiVir[p][sigma]*eta[lambda][nu];
3180+ }
3181+ value += twoElecInt*this->fockMatrix[moI][mu]*temp;
3182+ }
3183+ }
3184+ }
3185+ }
3186+ }
3187+ }
3188+ }
3189+ return value;
3190+}
3191+
3192+void ZindoS::CalcXiMatrices(double** xiOcc,
3193+ double** xiVir,
3194+ int exciteState,
3195+ double const* const* transposedFockMatrix) const{
3196+ int numberAOs = this->molecule->GetTotalNumberAOs();
3197+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3198+ int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
3199+ int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
3200+ MallocerFreer::GetInstance()->Initialize<double>(
3201+ xiOcc, numberActiveOcc, numberAOs);
3202+ MallocerFreer::GetInstance()->Initialize<double>(
3203+ xiVir, numberActiveVir, numberAOs);
3204+ stringstream ompErrors;
3205+ // xiOcc
3206+#pragma omp parallel for schedule(auto)
3207+ for(int p=0; p<numberActiveOcc; p++){
3208+ try{
3209+ for(int mu=0; mu<numberAOs; mu++){
3210+ for(int a=0; a<numberActiveVir; a++){
3211+ int moA = numberOcc + a;
3212+ int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(p,a);
3213+ xiOcc[p][mu] += this->matrixCIS[exciteState][slaterDeterminantIndex]
3214+ *transposedFockMatrix[mu][moA];
3215+ }
3216+ }
3217+ }
3218+ catch(MolDSException ex){
3219+#pragma omp critical
3220+ ompErrors << ex.what() << endl ;
3221+ }
3222+ }
3223+ // Exception throwing for omp-region
3224+ if(!ompErrors.str().empty()){
3225+ throw MolDSException(ompErrors.str());
3226+ }
3227+ // xiVir
3228+#pragma omp parallel for schedule(auto)
3229+ for(int p=0; p<numberActiveVir; p++){
3230+ try{
3231+ for(int mu=0; mu<numberAOs; mu++){
3232+ for(int i=0; i<numberActiveOcc; i++){
3233+ int moI = numberOcc - (i+1);
3234+ int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(i,p);
3235+ xiVir[p][mu] += this->matrixCIS[exciteState][slaterDeterminantIndex]
3236+ *transposedFockMatrix[mu][moI];
3237+ }
3238+ }
3239+ }
3240+ catch(MolDSException ex){
3241+#pragma omp critical
3242+ ompErrors << ex.what() << endl ;
3243+ }
3244+ }
3245+ // Exception throwing for omp-region
3246+ if(!ompErrors.str().empty()){
3247+ throw MolDSException(ompErrors.str());
3248+ }
3249+}
3250+
3251+// right hand side of (54) in [PT_1996]
3252+void ZindoS::CalcAuxiliaryVector(double* y,
3253+ double const* q,
3254+ double const* const* kRDagerGammaRInv,
3255+ const vector<MoIndexPair>& nonRedundantQIndeces,
3256+ const vector<MoIndexPair>& redundantQIndeces) const{
3257+ MallocerFreer::GetInstance()->Initialize<double>(
3258+ y,
3259+ nonRedundantQIndeces.size());
3260+ stringstream ompErrors;
3261+#pragma omp parallel for schedule(auto)
3262+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
3263+ try{
3264+ int moI = nonRedundantQIndeces[i].moI;
3265+ int moJ = nonRedundantQIndeces[i].moJ;
3266+ y[i] += q[i]/this->GetNNRElement(moI, moJ, moI, moJ);
3267+ for(int j=0; j<redundantQIndeces.size(); j++){
3268+ int k = nonRedundantQIndeces.size() + j;
3269+ y[i] += kRDagerGammaRInv[i][j]*q[k];
3270+ }
3271+ }
3272+ catch(MolDSException ex){
3273+#pragma omp critical
3274+ ompErrors << ex.what() << endl ;
3275+ }
3276+ }
3277+ // Exception throwing for omp-region
3278+ if(!ompErrors.str().empty()){
3279+ throw MolDSException(ompErrors.str());
3280+ }
3281+}
3282+
3283+// see (40) and (45) in [PT_1996].
3284+// This method calculates "\Gamma_{NR} - K_{NR}" to solve (54) in [PT_1966]
3285+// Note taht K_{NR} is not calculated.
3286+void ZindoS::CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR, const vector<MoIndexPair>& nonRedundantQIndeces) const{
3287+ stringstream ompErrors;
3288+#pragma omp parallel for schedule(auto)
3289+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
3290+ try{
3291+ int moI = nonRedundantQIndeces[i].moI;
3292+ int moJ = nonRedundantQIndeces[i].moJ;
3293+ for(int j=i; j<nonRedundantQIndeces.size(); j++){
3294+ int moK = nonRedundantQIndeces[j].moI;
3295+ int moL = nonRedundantQIndeces[j].moJ;
3296+ gammaNRMinusKNR[i][j] = this->GetGammaNRElement(moI, moJ, moK, moL)
3297+ -this->GetKNRElement(moI, moJ, moK, moL);
3298+ }
3299+ }
3300+ catch(MolDSException ex){
3301+#pragma omp critical
3302+ ompErrors << ex.what() << endl ;
3303+ }
3304+ }
3305+ // Exception throwing for omp-region
3306+ if(!ompErrors.str().empty()){
3307+ throw MolDSException(ompErrors.str());
3308+ }
3309+}
3310+
3311+// see (41), (42), and (46) in [PT_1996].
3312+// This method calculates "K_{R}^{\dagger} * Gamma_{R}" matrix, see (41), (42), and (46) to solve (54) in [PT_1996]
3313+// Note taht K_{R}^{\dager} is not calculated.
3314+void ZindoS::CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
3315+ const vector<MoIndexPair>& nonRedundantQIndeces,
3316+ const vector<MoIndexPair>& redundantQIndeces) const{
3317+ stringstream ompErrors;
3318+#pragma omp parallel for schedule(auto)
3319+ for(int i=0; i<nonRedundantQIndeces.size(); i++){
3320+ try{
3321+ int moI = nonRedundantQIndeces[i].moI;
3322+ int moJ = nonRedundantQIndeces[i].moJ;
3323+ for(int j=0; j<redundantQIndeces.size(); j++){
3324+ int moK = redundantQIndeces[j].moI;
3325+ int moL = redundantQIndeces[j].moJ;
3326+ kRDagerGammaRInv[i][j] = this->GetKRDagerElement(moI, moJ, moK, moL)
3327+ /this->GetGammaRElement(moK, moL, moK, moL);
3328+ }
3329+ }
3330+ catch(MolDSException ex){
3331+#pragma omp critical
3332+ ompErrors << ex.what() << endl ;
3333+ }
3334+ }
3335+ // Exception throwing for omp-region
3336+ if(!ompErrors.str().empty()){
3337+ throw MolDSException(ompErrors.str());
3338+ }
3339+}
3340+
3341+// see (40) in [PT_1996]
3342+double ZindoS::GetGammaNRElement(int moI, int moJ, int moK, int moL) const{
3343+ double value=0.0;
3344+ if(moI==moK && moJ==moL){
3345+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3346+ double nI = moI<numberOcc ? 2.0 : 0.0;
3347+ double nJ = moJ<numberOcc ? 2.0 : 0.0;
3348+ value = (this->energiesMO[moJ]-this->energiesMO[moI])/(nJ-nI);
3349+ }
3350+ return value;
3351+}
3352+
3353+// see (41) & (42) in [PT_1996]
3354+double ZindoS::GetGammaRElement(int moI, int moJ, int moK, int moL) const{
3355+ double value=0.0;
3356+ if(moI==moK && moJ==moL){
3357+ value = moI==moJ ? 1.0 : this->energiesMO[moJ]-this->energiesMO[moI];
3358+ }
3359+ return value;
3360+}
3361+
3362+// see (43) in [PT_1996]
3363+double ZindoS::GetNNRElement(int moI, int moJ, int moK, int moL) const{
3364+ double value=0.0;
3365+ if(moI==moK && moJ==moL){
3366+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3367+ double nI = moI<numberOcc ? 2.0 : 0.0;
3368+ double nJ = moJ<numberOcc ? 2.0 : 0.0;
3369+ value = (nJ-nI);
3370+ }
3371+ return value;
3372+}
3373+
3374+// see (44) in [PT_1996]
3375+double ZindoS::GetNRElement(int moI, int moJ, int moK, int moL) const{
3376+ double value=0.0;
3377+ if(moI==moK && moJ==moL){
3378+ value = 1.0;
3379+ }
3380+ return value;
3381+}
3382+
3383+// see (44) in [PT_1996]
3384+double ZindoS::GetKNRElement(int moI, int moJ, int moK, int moL) const{
3385+ double value=0.0;
3386+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3387+ int nI = moI<numberOcc ? 2 : 0;
3388+ int nJ = moJ<numberOcc ? 2 : 0;
3389+ int nK = moK<numberOcc ? 2 : 0;
3390+ int nL = moL<numberOcc ? 2 : 0;
3391+
3392+ if(nI!=nJ && nK!=nL){
3393+ value = this->GetAuxiliaryKNRKRElement(moI, moJ, moK, moL);
3394+ }
3395+ //See (24) in [DL_1990] about "0.5" multiplied to "GetKNRElement".
3396+ return 0.5*value;
3397+}
3398+
3399+// Dager of (45) in [PT_1996]. Note taht the (45) is real number.
3400+double ZindoS::GetKRDagerElement(int moI, int moJ, int moK, int moL) const{
3401+ return this->GetKRElement(moK, moL, moI, moJ);
3402+}
3403+
3404+// see (45) in [PT_1996]
3405+double ZindoS::GetKRElement(int moI, int moJ, int moK, int moL) const{
3406+ double value=0.0;
3407+ int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
3408+ int nI = moI<numberOcc ? 2 : 0;
3409+ int nJ = moJ<numberOcc ? 2 : 0;
3410+ int nK = moK<numberOcc ? 2 : 0;
3411+ int nL = moL<numberOcc ? 2 : 0;
3412+
3413+ if(nI==nJ && nK!=nL){
3414+ value = this->GetAuxiliaryKNRKRElement(moI, moJ, moK, moL);
3415+ }
3416+ //See (24) in [DL_1990] about "0.5" multiplied to "GetKRElement".
3417+ return 0.5*value;
3418+}
3419+
3420+double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
3421+ double value = 0.0;
3422+
3423+ // Fast algorith, but this is not easy to read.
3424+ // Slow algorithm is alos written below.
3425+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3426+ const Atom& atomA = *this->molecule->GetAtom(A);
3427+ int firstAOIndexA = atomA.GetFirstAOIndex();
3428+ int lastAOIndexA = atomA.GetLastAOIndex();
3429+
3430+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3431+ int muOffSet = mu - firstAOIndexA;
3432+ for(int nu=mu; nu<=lastAOIndexA; nu++){
3433+ int nuOffSet = nu - firstAOIndexA;
3434+ double tmpMN01 = 0.0, tmpMN02 = 0.0, tmpMN03 = 0.0,
3435+ tmpMN04 = 0.0, tmpMN05 = 0.0, tmpMN06 = 0.0,
3436+ tmpMN13 = 0.0, tmpMN14 = 0.0, tmpMN15 = 0.0,
3437+ tmpMN16 = 0.0, tmpMN17 = 0.0, tmpMN18 = 0.0;
3438+ tmpMN01 = 4.0
3439+ *this->fockMatrix[moI][mu]
3440+ *this->fockMatrix[moJ][nu];
3441+ tmpMN02 = 4.0
3442+ *this->fockMatrix[moK][mu]
3443+ *this->fockMatrix[moL][nu];
3444+ tmpMN03 = this->fockMatrix[moI][mu]
3445+ *this->fockMatrix[moK][nu];
3446+ tmpMN04 = this->fockMatrix[moJ][mu]
3447+ *this->fockMatrix[moL][nu];
3448+ tmpMN05 = this->fockMatrix[moI][mu]
3449+ *this->fockMatrix[moL][nu];
3450+ tmpMN06 = this->fockMatrix[moJ][mu]
3451+ *this->fockMatrix[moK][nu];
3452+ if(mu != nu){
3453+ tmpMN13 = 4.0
3454+ *this->fockMatrix[moI][nu]
3455+ *this->fockMatrix[moJ][mu];
3456+ tmpMN14 = 4.0
3457+ *this->fockMatrix[moK][nu]
3458+ *this->fockMatrix[moL][mu];
3459+ tmpMN15 = this->fockMatrix[moI][nu]
3460+ *this->fockMatrix[moK][mu];
3461+ tmpMN16 = this->fockMatrix[moJ][nu]
3462+ *this->fockMatrix[moL][mu];
3463+ tmpMN17 = this->fockMatrix[moI][nu]
3464+ *this->fockMatrix[moL][mu];
3465+ tmpMN18 = this->fockMatrix[moJ][nu]
3466+ *this->fockMatrix[moK][mu];
3467+ }
3468+
3469+ for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
3470+ const Atom& atomB = *this->molecule->GetAtom(B);
3471+ int firstAOIndexB = atomB.GetFirstAOIndex();
3472+ int lastAOIndexB = atomB.GetLastAOIndex();
3473+
3474+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3475+ int lambdaOffSet = lambda - firstAOIndexB;
3476+ double tmpMNL01 = 0.0, tmpMNL02 = 0.0, tmpMNL03 = 0.0, tmpMNL04 = 0.0,
3477+ tmpMNL05 = 0.0, tmpMNL06 = 0.0, tmpMNL07 = 0.0, tmpMNL08 = 0.0,
3478+ tmpMNL09 = 0.0, tmpMNL10 = 0.0, tmpMNL11 = 0.0, tmpMNL12 = 0.0,
3479+ tmpMNL13 = 0.0, tmpMNL14 = 0.0, tmpMNL15 = 0.0, tmpMNL16 = 0.0,
3480+ tmpMNL17 = 0.0, tmpMNL18 = 0.0, tmpMNL19 = 0.0, tmpMNL20 = 0.0,
3481+ tmpMNL21 = 0.0, tmpMNL22 = 0.0, tmpMNL23 = 0.0, tmpMNL24 = 0.0;
3482+ tmpMNL01 = tmpMN01*this->fockMatrix[moK][lambda];
3483+ tmpMNL02 = tmpMN02*this->fockMatrix[moI][lambda];
3484+ tmpMNL03 = tmpMN03*this->fockMatrix[moJ][lambda];
3485+ tmpMNL04 = tmpMN04*this->fockMatrix[moI][lambda];
3486+ tmpMNL05 = tmpMN05*this->fockMatrix[moJ][lambda];
3487+ tmpMNL06 = tmpMN06*this->fockMatrix[moI][lambda];
3488+ tmpMNL07 = tmpMN01*this->fockMatrix[moL][lambda];
3489+ tmpMNL08 = tmpMN02*this->fockMatrix[moJ][lambda];
3490+ tmpMNL09 = tmpMN03*this->fockMatrix[moL][lambda];
3491+ tmpMNL10 = tmpMN04*this->fockMatrix[moK][lambda];
3492+ tmpMNL11 = tmpMN05*this->fockMatrix[moK][lambda];
3493+ tmpMNL12 = tmpMN06*this->fockMatrix[moL][lambda];
3494+ tmpMNL01 -= tmpMNL03 + tmpMNL06;
3495+ tmpMNL04 += tmpMNL05;
3496+ tmpMNL08 -= tmpMNL10 + tmpMNL12;
3497+ tmpMNL09 += tmpMNL11;
3498+ if(mu != nu){
3499+ tmpMNL13 = tmpMN13*this->fockMatrix[moK][lambda];
3500+ tmpMNL14 = tmpMN14*this->fockMatrix[moI][lambda];
3501+ tmpMNL15 = tmpMN15*this->fockMatrix[moJ][lambda];
3502+ tmpMNL16 = tmpMN16*this->fockMatrix[moI][lambda];
3503+ tmpMNL17 = tmpMN17*this->fockMatrix[moJ][lambda];
3504+ tmpMNL18 = tmpMN18*this->fockMatrix[moI][lambda];
3505+ tmpMNL19 = tmpMN13*this->fockMatrix[moL][lambda];
3506+ tmpMNL20 = tmpMN14*this->fockMatrix[moJ][lambda];
3507+ tmpMNL21 = tmpMN15*this->fockMatrix[moL][lambda];
3508+ tmpMNL22 = tmpMN16*this->fockMatrix[moK][lambda];
3509+ tmpMNL23 = tmpMN17*this->fockMatrix[moK][lambda];
3510+ tmpMNL24 = tmpMN18*this->fockMatrix[moL][lambda];
3511+ tmpMNL13 -= tmpMNL15 + tmpMNL18;
3512+ tmpMNL16 += tmpMNL17;
3513+ tmpMNL20 -= tmpMNL22 + tmpMNL24;
3514+ tmpMNL21 += tmpMNL23;
3515+ tmpMNL01 += tmpMNL13;
3516+ tmpMNL02 += tmpMNL14;
3517+ tmpMNL04 += tmpMNL16;
3518+ tmpMNL07 += tmpMNL19;
3519+ tmpMNL08 += tmpMNL20;
3520+ tmpMNL09 += tmpMNL21;
3521+ }
3522+ for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
3523+ int sigmaOffSet = sigma - firstAOIndexB;
3524+ double tmpValue = 0.0;
3525+ tmpValue += tmpMNL01*this->fockMatrix[moL][sigma];
3526+ tmpValue += tmpMNL02*this->fockMatrix[moJ][sigma];
3527+ tmpValue -= tmpMNL04*this->fockMatrix[moK][sigma];
3528+ if(lambda != sigma){
3529+ tmpValue += tmpMNL07*this->fockMatrix[moK][sigma];
3530+ tmpValue += tmpMNL08*this->fockMatrix[moI][sigma];
3531+ tmpValue -= tmpMNL09*this->fockMatrix[moJ][sigma];
3532+ }
3533+ double gamma = 0.0;
3534+ if(A!=B){
3535+ gamma = this->twoElecTwoCore[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
3536+ }
3537+ else{
3538+ if(mu==nu && lambda==sigma){
3539+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3540+ OrbitalType orbitalLambda = atomA.GetValence(lambdaOffSet);
3541+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3542+ }
3543+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3544+ OrbitalType orbitalMu = atomA.GetValence(muOffSet);
3545+ OrbitalType orbitalNu = atomA.GetValence(nuOffSet);
3546+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3547+ }
3548+ else{
3549+ gamma = 0.0;
3550+ }
3551+ gamma *= 0.5;
3552+ }
3553+ value += tmpValue*gamma;
3554+ }
3555+ }
3556+ }
3557+ }
3558+ }
3559+ }
3560+ // End of the fast algorith.
3561+
3562+ /*
3563+ // Algorithm using blas
3564+ double** twoElec = NULL;
3565+ double* twiceMoIJ = NULL;
3566+ double* twiceMoIK = NULL;
3567+ double* twiceMoIL = NULL;
3568+ double* twiceMoKL = NULL;
3569+ double* twiceMoJL = NULL;
3570+ double* twiceMoJK = NULL;
3571+ double* tmpVector = NULL;
3572+ int numAOs = this->molecule->GetTotalNumberAOs();
3573+ MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3574+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3575+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3576+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3577+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
3578+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
3579+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
3580+ MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
3581+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3582+ const Atom& atomA = *this->molecule->GetAtom(A);
3583+ int firstAOIndexA = atomA.GetFirstAOIndex();
3584+ int lastAOIndexA = atomA.GetLastAOIndex();
3585+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3586+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3587+ twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
3588+ twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
3589+ twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
3590+ }
3591+ }
3592+ }
3593+
3594+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3595+ const Atom& atomB = *this->molecule->GetAtom(B);
3596+ int firstAOIndexB = atomB.GetFirstAOIndex();
3597+ int lastAOIndexB = atomB.GetLastAOIndex();
3598+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3599+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3600+ twiceMoKL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
3601+ twiceMoJL[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
3602+ twiceMoJK[B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
3603+ }
3604+ }
3605+ }
3606+
3607+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3608+ const Atom& atomA = *this->molecule->GetAtom(A);
3609+ int firstAOIndexA = atomA.GetFirstAOIndex();
3610+ int lastAOIndexA = atomA.GetLastAOIndex();
3611+ for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
3612+ const Atom& atomB = *this->molecule->GetAtom(B);
3613+ int firstAOIndexB = atomB.GetFirstAOIndex();
3614+ int lastAOIndexB = atomB.GetLastAOIndex();
3615+ double gamma = 0.0;
3616+ if(A!=B){
3617+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3618+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3619+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3620+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3621+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3622+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3623+ this->twoElecTwoCore[A]
3624+ [B]
3625+ [mu-firstAOIndexA]
3626+ [nu-firstAOIndexA]
3627+ [lambda-firstAOIndexB]
3628+ [sigma-firstAOIndexB];
3629+ }
3630+ }
3631+ }
3632+ }
3633+ }
3634+ else{
3635+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3636+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3637+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3638+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3639+ if(mu==nu && lambda==sigma){
3640+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3641+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3642+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3643+ }
3644+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3645+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3646+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
3647+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3648+ }
3649+ else{
3650+ gamma = 0.0;
3651+ }
3652+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3653+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
3654+ }
3655+ }
3656+ }
3657+ }
3658+ }
3659+ }
3660+ }
3661+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3662+ twoElec,
3663+ twiceMoKL,
3664+ tmpVector);
3665+ value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, tmpVector);
3666+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3667+ twoElec,
3668+ twiceMoJL,
3669+ tmpVector);
3670+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, tmpVector);
3671+ MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy,
3672+ twoElec,
3673+ twiceMoJK,
3674+ tmpVector);
3675+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, tmpVector);
3676+ MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3677+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3678+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3679+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3680+ MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
3681+ MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
3682+ MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
3683+ MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
3684+ // End of algorithm using blas
3685+ */
3686+
3687+ /*
3688+ // Second algorithm using blas.
3689+ // This algorithm uses DGEMM.
3690+ double** twoElec = NULL;
3691+ double* twiceMoIJ = NULL;
3692+ double* twiceMoIK = NULL;
3693+ double* twiceMoIL = NULL;
3694+ double** twiceMoB = NULL;
3695+ double** tmpMatrix = NULL;
3696+ int numAOs = this->molecule->GetTotalNumberAOs();
3697+ MallocerFreer::GetInstance()->Malloc<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3698+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3699+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3700+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3701+ MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
3702+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
3703+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3704+ const Atom& atomA = *this->molecule->GetAtom(A);
3705+ int firstAOIndexA = atomA.GetFirstAOIndex();
3706+ int lastAOIndexA = atomA.GetLastAOIndex();
3707+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3708+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3709+ twiceMoIJ[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moJ][nu ];
3710+ twiceMoIK[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moK][nu ];
3711+ twiceMoIL[A*dxy*dxy+(mu -firstAOIndexA)*dxy+(nu -firstAOIndexA)]=fockMatrix[moI][mu ]*fockMatrix[moL][nu ];
3712+ }
3713+ }
3714+ }
3715+
3716+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3717+ const Atom& atomB = *this->molecule->GetAtom(B);
3718+ int firstAOIndexB = atomB.GetFirstAOIndex();
3719+ int lastAOIndexB = atomB.GetLastAOIndex();
3720+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3721+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3722+ twiceMoB[0][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moK][lambda]*fockMatrix[moL][sigma];
3723+ twiceMoB[1][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moL][sigma];
3724+ twiceMoB[2][B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)]=fockMatrix[moJ][lambda]*fockMatrix[moK][sigma];
3725+ }
3726+ }
3727+ }
3728+
3729+ for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
3730+ const Atom& atomA = *this->molecule->GetAtom(A);
3731+ int firstAOIndexA = atomA.GetFirstAOIndex();
3732+ int lastAOIndexA = atomA.GetLastAOIndex();
3733+ for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
3734+ const Atom& atomB = *this->molecule->GetAtom(B);
3735+ int firstAOIndexB = atomB.GetFirstAOIndex();
3736+ int lastAOIndexB = atomB.GetLastAOIndex();
3737+ double gamma = 0.0;
3738+ if(A!=B){
3739+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3740+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3741+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3742+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3743+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3744+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3745+ this->twoElecTwoCore[A]
3746+ [B]
3747+ [mu-firstAOIndexA]
3748+ [nu-firstAOIndexA]
3749+ [lambda-firstAOIndexB]
3750+ [sigma-firstAOIndexB];
3751+ }
3752+ }
3753+ }
3754+ }
3755+ }
3756+ else{
3757+ for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
3758+ for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
3759+ for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
3760+ for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
3761+ if(mu==nu && lambda==sigma){
3762+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3763+ OrbitalType orbitalLambda = atomB.GetValence(lambda-firstAOIndexB);
3764+ gamma = this->GetCoulombInt(orbitalMu, orbitalLambda, atomA);
3765+ }
3766+ else if((mu==lambda && nu==sigma) || (nu==lambda && mu==sigma) ){
3767+ OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
3768+ OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
3769+ gamma = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
3770+ }
3771+ else{
3772+ gamma = 0.0;
3773+ }
3774+ twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
3775+ [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] = gamma;
3776+ }
3777+ }
3778+ }
3779+ }
3780+ }
3781+ }
3782+ }
3783+
3784+ MolDS_wrappers::Blas::GetInstance()->Dgemm(false, true, true,
3785+ this->molecule->GetNumberAtoms()*dxy*dxy,
3786+ 3,
3787+ this->molecule->GetNumberAtoms()*dxy*dxy,
3788+ 1.0,
3789+ twoElec,
3790+ twiceMoB,
3791+ 0.0,
3792+ tmpMatrix);
3793+ value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
3794+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
3795+ value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
3796+ MallocerFreer::GetInstance()->Free<double>(&twoElec, this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
3797+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
3798+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
3799+ MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
3800+ MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
3801+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
3802+ // End of second algorithm using blas
3803+ */
3804+
3805+ /*
3806+ // slow algorithm
3807+ value = 4.0*this->GetMolecularIntegralElement(moI, moJ, moK, moL,
3808+ *this->molecule,
3809+ this->fockMatrix, NULL)
3810+ -1.0*this->GetMolecularIntegralElement(moI, moK, moJ, moL,
3811+ *this->molecule,
3812+ this->fockMatrix, NULL)
3813+ -1.0*this->GetMolecularIntegralElement(moI, moL, moJ, moK,
3814+ *this->molecule,
3815+ this->fockMatrix, NULL);
3816+ */
3817+ return value;
3818+}
3819+
25503820 // elecStates is indeces of the electroinc eigen states.
25513821 // The index = 0 means electronic ground state.
25523822 void ZindoS::CalcForce(const vector<int>& elecStates){
@@ -2734,5 +4004,3 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
27344004
27354005 }
27364006
2737-
2738-
--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -45,6 +45,84 @@ protected:
4545 std::string messageOmpElapsedTimeCalcCISMarix;
4646 std::string messageOmpElapsedTimeCIS;
4747 std::string messageDoneCalcCISMatrix;
48+ /*** from MNDO ***/
49+ std::string errorMessageCalcZMatrixForceEtaNull;
50+ int zMatrixForceElecStatesNum;
51+ int etaMatrixForceElecStatesNum;
52+ double*** zMatrixForce;
53+ double*** etaMatrixForce;
54+ struct MoIndexPair{int moI; int moJ; bool isMoICIMO; bool isMoJCIMO;};
55+ void CheckZMatrixForce(const std::vector<int>& elecStates);
56+ void CheckEtaMatrixForce(const std::vector<int>& elecStates);
57+ void CalcEtaMatrixForce(const std::vector<int>& elecStates);
58+ void CalcZMatrixForce(const std::vector<int>& elecStates);
59+ double GetZMatrixForceElement(double const* y,
60+ double const* q,
61+ double const* const* transposedFockMatrix,
62+ const std::vector<MoIndexPair>& nonRedundantQIndeces,
63+ const std::vector<MoIndexPair>& redundantQIndeces,
64+ int mu,
65+ int nu) const;
66+ void MallocTempMatrixForZMatrix(double** delta,
67+ double** q,
68+ double*** gammaNRMinusKNR,
69+ double*** kRDag,
70+ double** y,
71+ double*** transposedFockMatrix,
72+ double*** xiOcc,
73+ double*** xiVir,
74+ int sizeQNR,
75+ int sizeQR) const;
76+ void FreeTempMatrixForZMatrix(double** delta,
77+ double** q,
78+ double*** gammaNRMinusKNR,
79+ double*** kRDag,
80+ double** y,
81+ double*** transposedFockMatrix,
82+ double*** xiOcc,
83+ double*** xiVir,
84+ int sizeQNR,
85+ int sizeQR) const;
86+ void CalcDeltaVector(double* delta, int exciteState) const;
87+ void CalcActiveSetVariablesQ(std::vector<MoIndexPair>* nonRedundantQIndeces,
88+ std::vector<MoIndexPair>* redundantQIndeces,
89+ int numberActiveOcc,
90+ int numberActiveVir) const;
91+ void CalcQVector(double* q,
92+ double const* delta,
93+ double const* const* xiOcc,
94+ double const* const* xiVir,
95+ double const* const* eta,
96+ const std::vector<MoIndexPair>& nonRedundantQIndeces,
97+ 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;
103+ void CalcXiMatrices(double** xiOcc,
104+ double** xiVir,
105+ int exciteState,
106+ double const* const* transposedFockMatrix) const;
107+ void CalcAuxiliaryVector(double* y,
108+ double const* q,
109+ double const* const* kRDagerGammaRInv,
110+ const std::vector<MoIndexPair>& nonRedundantQIndeces,
111+ const std::vector<MoIndexPair>& redundantQIndeces) const;
112+ double GetGammaNRElement(int moI, int moJ, int moK, int moL) const;
113+ double GetGammaRElement(int moI, int moJ, int moK, int moL) const;
114+ double GetNNRElement(int moI, int moJ, int moK, int moL) const;
115+ double GetNRElement(int moI, int moJ, int moK, int moL) const;
116+ double GetKNRElement(int moI, int moJ, int moK, int moL) const;
117+ double GetKRElement(int moI, int moJ, int moK, int moL) const;
118+ double GetKRDagerElement(int moI, int moJ, int moK, int moL) const;
119+ double GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const;
120+ void CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR,
121+ const std::vector<MoIndexPair>& nonRedundantQIndeces) const;
122+ void CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
123+ const std::vector<MoIndexPair>& nonRedundantQIndeces,
124+ const std::vector<MoIndexPair>& redundantQIndeces) const;
125+ /*** end from MNDO ***/
48126 virtual void SetMessages();
49127 virtual void SetEnableAtomTypes();
50128 virtual void CalcCISProperties();
@@ -114,6 +192,7 @@ protected:
114192 int moA,
115193 int moJ,
116194 int moB) const;
195+ bool RequiresExcitedStatesForce(const std::vector<int>& elecStates) const;
117196 virtual void CalcForce(const std::vector<int>& elecStates);
118197 int GetSlaterDeterminantIndex(int activeOccIndex, int activeVirIndex) const;
119198 int GetActiveOccIndex(const MolDS_base::Molecule& molecule, int matrixCISIndex) const;