OSDN -- Open-Source-Software-Entwicklung und Downloads
Übersetzungsstatus von Deutsch translation status for de
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
Keine Tags

Frequently used words (click to add to your profile)

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

Commit MetaInfo

Revision4c222a5c90c2b3cd72c9585170a8b1973f3e9cd8 (tree)
Zeit2011-07-16 20:49:31
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Speed up of ZindoS::CalcCISMatrix is carried out.

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

Ändern Zusammenfassung

Diff

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -75,6 +75,7 @@ private:
7575 string errorMessageDavidsonNotConverged;
7676 string errorMessageDavidsonMaxIter;
7777 string errorMessageDavidsonMaxDim;
78+ string errorMessageCalcCISMatrix;
7879 string messageStartCalcCISMatrix;
7980 string messageConsumedCalcCISMatrix;
8081 string messageDoneCalcCISMatrix;
@@ -129,6 +130,8 @@ void ZindoS::SetMessages(){
129130 this->errorMessageNishimotoMataga = "Error in base_zindo::ZindoS::GetNishimotoMatagaTwoEleInt: Invalid orbitalType.\n";
130131 this->errorMessageMolecularIntegralElement
131132 = "Error in zindo::ZindoS::GetMolecularIntegralElement: Non available orbital is contained.\n";
133+ this->errorMessageCalcCISMatrix
134+ = "Error in zindo::ZindoS::CalcCISMatrix: Non available orbital is contained.\n";
132135 this->errorMessageDavidsonNotConverged = "Error in zindo::ZindoS::DoesCISDavidson: Davidson did not met convergence criterion. \n";
133136 this->errorMessageDavidsonMaxIter = "Davidson roop reaches max_iter=";
134137 this->errorMessageDavidsonMaxDim = "Dimension of the expansion vectors reaches max_dim=";
@@ -618,7 +621,7 @@ double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
618621 orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
619622
620623 // CNDO term
621- for(int B=0; B<molecule->GetAtomVect()->size(); B++){
624+ for(int B=A; B<molecule->GetAtomVect()->size(); B++){
622625 atomB = (*molecule->GetAtomVect())[B];
623626 firstAOIndexB = atomB->GetFirstAOIndex();
624627 numberAOsB = atomB->GetValence().size();
@@ -626,14 +629,16 @@ double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
626629 for(int nu=firstAOIndexB; nu<firstAOIndexB+numberAOsB; nu++){
627630 orbitalNu = atomB->GetValence()[nu-firstAOIndexB];
628631
629- if(A==B){
630- gamma = atomA->GetZindoF0ss();
632+ if(A<B){
633+ gamma = this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalNu);
634+ value += gamma*fockMatrix[moI][mu]*fockMatrix[moJ][mu]*fockMatrix[moK][nu]*fockMatrix[moL][nu];
635+ value += gamma*fockMatrix[moI][nu]*fockMatrix[moJ][nu]*fockMatrix[moK][mu]*fockMatrix[moL][mu];
631636 }
632637 else{
633- gamma = this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalNu);
638+ gamma = atomA->GetZindoF0ss();
639+ value += gamma*fockMatrix[moI][mu]*fockMatrix[moJ][mu]*fockMatrix[moK][nu]*fockMatrix[moL][nu];
634640 }
635641
636- value += gamma*fockMatrix[moI][mu]*fockMatrix[moJ][mu]*fockMatrix[moK][nu]*fockMatrix[moL][nu];
637642 }
638643 }
639644
@@ -1044,12 +1049,235 @@ void ZindoS::CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir){
10441049 int moI = this->molecule->GetTotalNumberValenceElectrons()/2 - (k/numberVir) -1;
10451050 int moA = this->molecule->GetTotalNumberValenceElectrons()/2 + (k%numberVir);
10461051
1047- for(int l=0; l<numberOcc*numberVir; l++){
1052+ for(int l=k; l<numberOcc*numberVir; l++){
10481053 // single excitation from J-th (occupied)MO to B-th (virtual)MO
10491054 int moJ = this->molecule->GetTotalNumberValenceElectrons()/2 - (l/numberVir) -1;
10501055 int moB = this->molecule->GetTotalNumberValenceElectrons()/2 + (l%numberVir);
1051-
10521056 double value=0.0;
1057+
1058+ // Fast algorith, but this is not easy to read. Slow algorithm is alos written below.
1059+ Atom* atomA = NULL;
1060+ Atom* atomB = NULL;
1061+ int firstAOIndexA;
1062+ int firstAOIndexB;
1063+ int numberAOsA;
1064+ int numberAOsB;
1065+ double gamma;
1066+ double exchange;
1067+ double coulomb;
1068+ OrbitalType orbitalMu;
1069+ OrbitalType orbitalNu;
1070+ // Off diagonal term (right upper)
1071+ if(k<l){
1072+ for(int A=0; A<molecule->GetAtomVect()->size(); A++){
1073+ atomA = (*molecule->GetAtomVect())[A];
1074+ firstAOIndexA = atomA->GetFirstAOIndex();
1075+ numberAOsA = atomA->GetValence().size();
1076+
1077+ for(int mu=firstAOIndexA; mu<firstAOIndexA+numberAOsA; mu++){
1078+ orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
1079+
1080+ // CNDO term
1081+ for(int B=A; B<molecule->GetAtomVect()->size(); B++){
1082+ atomB = (*molecule->GetAtomVect())[B];
1083+ firstAOIndexB = atomB->GetFirstAOIndex();
1084+ numberAOsB = atomB->GetValence().size();
1085+
1086+ for(int nu=firstAOIndexB; nu<firstAOIndexB+numberAOsB; nu++){
1087+ orbitalNu = atomB->GetValence()[nu-firstAOIndexB];
1088+
1089+ if(A<B){
1090+ gamma = this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalNu);
1091+ value += 2.0*gamma*fockMatrix[moA][mu]
1092+ *fockMatrix[moI][mu]
1093+ *fockMatrix[moJ][nu]
1094+ *fockMatrix[moB][nu];
1095+ value -= gamma*fockMatrix[moA][mu]
1096+ *fockMatrix[moB][mu]
1097+ *fockMatrix[moI][nu]
1098+ *fockMatrix[moJ][nu];
1099+ value += 2.0*gamma*fockMatrix[moA][nu]
1100+ *fockMatrix[moI][nu]
1101+ *fockMatrix[moJ][mu]
1102+ *fockMatrix[moB][mu];
1103+ value -= gamma*fockMatrix[moA][nu]
1104+ *fockMatrix[moB][nu]
1105+ *fockMatrix[moI][mu]
1106+ *fockMatrix[moJ][mu];
1107+ }
1108+ else{
1109+ gamma = atomA->GetZindoF0ss();
1110+ value += 2.0*gamma*fockMatrix[moA][mu]
1111+ *fockMatrix[moI][mu]
1112+ *fockMatrix[moJ][nu]
1113+ *fockMatrix[moB][nu];
1114+ value -= gamma*fockMatrix[moA][mu]
1115+ *fockMatrix[moB][mu]
1116+ *fockMatrix[moI][nu]
1117+ *fockMatrix[moJ][nu];
1118+ }
1119+ }
1120+ }
1121+
1122+ // Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
1123+ for(int nu=firstAOIndexA; nu<firstAOIndexA+numberAOsA; nu++){
1124+ orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
1125+
1126+ if(mu!=nu){
1127+ exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
1128+ value += 2.0*exchange*fockMatrix[moA][mu]
1129+ *fockMatrix[moI][nu]
1130+ *fockMatrix[moJ][nu]
1131+ *fockMatrix[moB][mu];
1132+ value += 2.0*exchange*fockMatrix[moA][mu]
1133+ *fockMatrix[moI][nu]
1134+ *fockMatrix[moJ][mu]
1135+ *fockMatrix[moB][nu];
1136+ value -= exchange*fockMatrix[moA][mu]
1137+ *fockMatrix[moB][nu]
1138+ *fockMatrix[moI][nu]
1139+ *fockMatrix[moJ][mu];
1140+ value -= exchange*fockMatrix[moA][mu]
1141+ *fockMatrix[moB][nu]
1142+ *fockMatrix[moI][mu]
1143+ *fockMatrix[moJ][nu];
1144+ }
1145+
1146+ coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
1147+
1148+ if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
1149+ (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
1150+ gamma = atomA->GetZindoF0ss();
1151+ }
1152+ else{
1153+ stringstream ss;
1154+ ss << this->errorMessageCalcCISMatrix;
1155+ ss << this->errorMessageAtomType << AtomTypeStr(atomA->GetAtomType()) << "\n";
1156+ ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
1157+ ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
1158+ throw MolDSException(ss.str());
1159+ }
1160+
1161+ value += 2.0*(coulomb-gamma)*fockMatrix[moA][mu]
1162+ *fockMatrix[moI][mu]
1163+ *fockMatrix[moJ][nu]
1164+ *fockMatrix[moB][nu];
1165+ value -= (coulomb-gamma)*fockMatrix[moA][mu]
1166+ *fockMatrix[moB][mu]
1167+ *fockMatrix[moI][nu]
1168+ *fockMatrix[moJ][nu];
1169+ }
1170+ }
1171+ }
1172+ }
1173+ // Diagonal term
1174+ else if(k==l){
1175+ value = this->energiesMO[moA] - this->energiesMO[moI];
1176+ for(int A=0; A<molecule->GetAtomVect()->size(); A++){
1177+ atomA = (*molecule->GetAtomVect())[A];
1178+ firstAOIndexA = atomA->GetFirstAOIndex();
1179+ numberAOsA = atomA->GetValence().size();
1180+
1181+ for(int mu=firstAOIndexA; mu<firstAOIndexA+numberAOsA; mu++){
1182+ orbitalMu = atomA->GetValence()[mu-firstAOIndexA];
1183+
1184+ // CNDO term
1185+ for(int B=A; B<molecule->GetAtomVect()->size(); B++){
1186+ atomB = (*molecule->GetAtomVect())[B];
1187+ firstAOIndexB = atomB->GetFirstAOIndex();
1188+ numberAOsB = atomB->GetValence().size();
1189+
1190+ for(int nu=firstAOIndexB; nu<firstAOIndexB+numberAOsB; nu++){
1191+ orbitalNu = atomB->GetValence()[nu-firstAOIndexB];
1192+
1193+ if(A<B){
1194+ gamma = this->GetNishimotoMatagaTwoEleInt(atomA, orbitalMu, atomB, orbitalNu);
1195+ value += 2.0*gamma*fockMatrix[moI][mu]
1196+ *fockMatrix[moA][mu]
1197+ *fockMatrix[moA][nu]
1198+ *fockMatrix[moI][nu];
1199+ value -= gamma*fockMatrix[moI][mu]
1200+ *fockMatrix[moI][mu]
1201+ *fockMatrix[moA][nu]
1202+ *fockMatrix[moA][nu];
1203+ value += 2.0*gamma*fockMatrix[moI][nu]
1204+ *fockMatrix[moA][nu]
1205+ *fockMatrix[moA][mu]
1206+ *fockMatrix[moI][mu];
1207+ value -= gamma*fockMatrix[moI][nu]
1208+ *fockMatrix[moI][nu]
1209+ *fockMatrix[moA][mu]
1210+ *fockMatrix[moA][mu];
1211+ }
1212+ else{
1213+ gamma = atomA->GetZindoF0ss();
1214+ value += 2.0*gamma*fockMatrix[moI][mu]
1215+ *fockMatrix[moA][mu]
1216+ *fockMatrix[moA][nu]
1217+ *fockMatrix[moI][nu];
1218+ value -= gamma*fockMatrix[moI][mu]
1219+ *fockMatrix[moI][mu]
1220+ *fockMatrix[moA][nu]
1221+ *fockMatrix[moA][nu];
1222+ }
1223+ }
1224+ }
1225+
1226+ // Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
1227+ for(int nu=firstAOIndexA; nu<firstAOIndexA+numberAOsA; nu++){
1228+ orbitalNu = atomA->GetValence()[nu-firstAOIndexA];
1229+
1230+ if(mu!=nu){
1231+ exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
1232+ value += 2.0*exchange*fockMatrix[moI][mu]
1233+ *fockMatrix[moA][nu]
1234+ *fockMatrix[moA][nu]
1235+ *fockMatrix[moI][mu];
1236+ value += 2.0*exchange*fockMatrix[moI][mu]
1237+ *fockMatrix[moA][nu]
1238+ *fockMatrix[moA][mu]
1239+ *fockMatrix[moI][nu];
1240+ value -= exchange*fockMatrix[moI][mu]
1241+ *fockMatrix[moI][nu]
1242+ *fockMatrix[moA][nu]
1243+ *fockMatrix[moA][mu];
1244+ value -= exchange*fockMatrix[moI][mu]
1245+ *fockMatrix[moI][nu]
1246+ *fockMatrix[moA][mu]
1247+ *fockMatrix[moA][nu];
1248+ }
1249+
1250+ coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
1251+
1252+ if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
1253+ (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
1254+ gamma = atomA->GetZindoF0ss();
1255+ }
1256+ else{
1257+ stringstream ss;
1258+ ss << this->errorMessageCalcCISMatrix;
1259+ ss << this->errorMessageAtomType << AtomTypeStr(atomA->GetAtomType()) << "\n";
1260+ ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
1261+ ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
1262+ throw MolDSException(ss.str());
1263+ }
1264+
1265+ value += 2.0*(coulomb-gamma)*fockMatrix[moI][mu]
1266+ *fockMatrix[moA][mu]
1267+ *fockMatrix[moA][nu]
1268+ *fockMatrix[moI][nu];
1269+ value -= (coulomb-gamma)*fockMatrix[moI][mu]
1270+ *fockMatrix[moI][mu]
1271+ *fockMatrix[moA][nu]
1272+ *fockMatrix[moA][nu];
1273+ }
1274+ }
1275+ }
1276+ }
1277+ // End of the fast algorith.
1278+
1279+ /*
1280+ // Slow algorith, but this is easy to read. Fast altorithm is also written above.
10531281 // diagonal term
10541282 if(k==l){
10551283 value = this->energiesMO[moA] - this->energiesMO[moI]
@@ -1059,14 +1287,16 @@ void ZindoS::CalcCISMatrix(double** matrixCIS, int numberOcc, int numberVir){
10591287 this->molecule, this->fockMatrix, NULL);
10601288
10611289 }
1062- // off diagonal term (right upper)
1290+ // Off diagonal term (right upper)
10631291 else if(k<l){
10641292 value = 2.0*this->GetMolecularIntegralElement(moA, moI, moJ, moB,
10651293 this->molecule, this->fockMatrix, NULL)
10661294 - this->GetMolecularIntegralElement(moA, moB, moI, moJ,
10671295 this->molecule, this->fockMatrix, NULL);
10681296 }
1069-
1297+ // End of the slow algorith.
1298+ */
1299+
10701300 matrixCIS[k][l] = value;
10711301 }
10721302 }