修订版 | 10a3740bcf663a52b5a1ccabdd39c8bbac28d607 (tree) |
---|---|
时间 | 2012-12-03 18:45:00 |
作者 | Mikiya Fujii <mikiya.fujii@gmai...> |
Commiter | Mikiya Fujii |
mulliken population for excited states. #28605
git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1165 1136aad2-a195-0410-b898-f5ea1d11b9d8
@@ -134,6 +134,7 @@ void InputParser::SetMessages(){ | ||
134 | 134 | this->messageCisExcitonEnergies = "\t\tExciton energies: "; |
135 | 135 | this->messageCisAllTransitionDipoleMoments = "\t\tAll transition dipole moments: "; |
136 | 136 | this->messageCisNumPrintCoefficients = "\t\tNumber of printed coefficients of CIS-eigenvector: "; |
137 | + this->messageCisMulliken = "\t\tMulliken population of excited states: "; | |
137 | 138 | |
138 | 139 | // memory |
139 | 140 | this->messageMemoryConditions = "\tMemory conditions:\n"; |
@@ -309,6 +310,7 @@ void InputParser::SetMessages(){ | ||
309 | 310 | this->stringCISExcitonEnergies = "exciton_energies"; |
310 | 311 | this->stringCISAllTransitionDipoleMoments = "all_transition_dipole_moments"; |
311 | 312 | this->stringCISNumPrintCoefficients = "num_print_coefficients"; |
313 | + this->stringCISMulliken = "mulliken"; | |
312 | 314 | |
313 | 315 | // Memory |
314 | 316 | this->stringMemory = "memory"; |
@@ -812,6 +814,15 @@ int InputParser::ParseConditionsCIS(vector<string>* inputTerms, int parseIndex) | ||
812 | 814 | } |
813 | 815 | parseIndex++; |
814 | 816 | } |
817 | + // mulliken | |
818 | + if((*inputTerms)[parseIndex].compare(this->stringCISMulliken) == 0){ | |
819 | + int groundStateIndex = 0; | |
820 | + int elecStateIndex = atoi((*inputTerms)[parseIndex+1].c_str()); | |
821 | + if(groundStateIndex<elecStateIndex){ | |
822 | + Parameters::GetInstance()->AddElectronicStateIndexMullikenCIS(elecStateIndex); | |
823 | + } | |
824 | + parseIndex++; | |
825 | + } | |
815 | 826 | parseIndex++; |
816 | 827 | } |
817 | 828 | return parseIndex; |
@@ -1326,6 +1337,20 @@ void InputParser::ValidateCisConditions(const Molecule& molecule) const{ | ||
1326 | 1337 | Parameters::GetInstance()->SetNumberPrintCoefficientsCIS(numberSlaterDeterminants); |
1327 | 1338 | } |
1328 | 1339 | |
1340 | + // Validate electronic state for Mulliken population analysis | |
1341 | + if(Parameters::GetInstance()->RequiresMullikenCIS()){ | |
1342 | + vector<int>* indecesMulliken = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS(); | |
1343 | + int numExcitedStates = Parameters::GetInstance()->GetNumberExcitedStatesCIS(); | |
1344 | + vector<int>::iterator it=(*indecesMulliken).begin(); | |
1345 | + vector<int>::iterator end=(*indecesMulliken).end(); | |
1346 | + while(it<end){ | |
1347 | + if(numExcitedStates<*it){ | |
1348 | + it = (*indecesMulliken).erase(it); | |
1349 | + end = (*indecesMulliken).end(); | |
1350 | + } | |
1351 | + ++it; | |
1352 | + } | |
1353 | + } | |
1329 | 1354 | } |
1330 | 1355 | |
1331 | 1356 | void InputParser::ValidateMdConditions(const Molecule& molecule) const{ |
@@ -1585,6 +1610,15 @@ void InputParser::OutputCisConditions() const{ | ||
1585 | 1610 | } |
1586 | 1611 | this->OutputLog("\n"); |
1587 | 1612 | |
1613 | + if(Parameters::GetInstance()->RequiresMullikenCIS()){ | |
1614 | + vector<int>* indeces = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS(); | |
1615 | + for(int i=0; i<indeces->size(); i++){ | |
1616 | + this->OutputLog(boost::format("%s%d\n") % this->messageCisMulliken.c_str() | |
1617 | + % (*indeces)[i]); | |
1618 | + } | |
1619 | + } | |
1620 | + this->OutputLog("\n"); | |
1621 | + | |
1588 | 1622 | this->OutputLog("\n"); |
1589 | 1623 | } |
1590 | 1624 |
@@ -80,6 +80,7 @@ private: | ||
80 | 80 | std::string messageCisExcitonEnergies; |
81 | 81 | std::string messageCisAllTransitionDipoleMoments; |
82 | 82 | std::string messageCisNumPrintCoefficients; |
83 | + std::string messageCisMulliken; | |
83 | 84 | // Memory |
84 | 85 | std::string messageMemoryConditions; |
85 | 86 | std::string messageMemoryLimitHeap; |
@@ -234,6 +235,7 @@ private: | ||
234 | 235 | std::string stringCISExcitonEnergies; |
235 | 236 | std::string stringCISAllTransitionDipoleMoments; |
236 | 237 | std::string stringCISNumPrintCoefficients; |
238 | + std::string stringCISMulliken; | |
237 | 239 | // Memory |
238 | 240 | std::string stringMemory; |
239 | 241 | std::string stringMemoryEnd; |
@@ -61,6 +61,7 @@ Parameters::Parameters(){ | ||
61 | 61 | this->indecesMOPlot = NULL; |
62 | 62 | this->elecIndecesHolePlot = NULL; |
63 | 63 | this->elecIndecesParticlePlot = NULL; |
64 | + this->electronicStateIndecesMullikenCIS = NULL; | |
64 | 65 | } |
65 | 66 | |
66 | 67 | Parameters::~Parameters(){ |
@@ -83,6 +84,11 @@ Parameters::~Parameters(){ | ||
83 | 84 | this->elecIndecesParticlePlot = NULL; |
84 | 85 | //this->OutputLog("elecIndecesParticlePlot deleted\n"); |
85 | 86 | } |
87 | + if(this->electronicStateIndecesMullikenCIS != NULL){ | |
88 | + delete this->electronicStateIndecesMullikenCIS; | |
89 | + this->electronicStateIndecesMullikenCIS= NULL; | |
90 | + this->OutputLog("electronicStateIndecesMullikenCIS deleted\n"); | |
91 | + } | |
86 | 92 | } |
87 | 93 | |
88 | 94 | Parameters* Parameters::GetInstance(){ |
@@ -211,6 +217,8 @@ void Parameters::SetMessages(){ | ||
211 | 217 | = "Error in base::Parameters::GetIndecesHolePlot: elecIndecesHolePlot is NULL.\n"; |
212 | 218 | this->errorMessageGetIndecesParticlePlotNull |
213 | 219 | = "Error in base::Parameters::GetIndecesParticlePlot: elecIndecesParticlePlot is NULL.\n"; |
220 | + this->errorMessageGetElectronicStateIndecesMullikenCISNull | |
221 | + = "Error in base::Parameters::GetElectronicStateIndecesMullikenCIS: electronicStateIndecesMullikenCIS is NULL.\n"; | |
214 | 222 | } |
215 | 223 | |
216 | 224 | SimulationType Parameters::GetCurrentSimulation() const{ |
@@ -675,6 +683,27 @@ void Parameters::SetRequiresAllTransitionDipoleMomentsCIS(bool requiresAllTransi | ||
675 | 683 | this->requiresAllTransitionDipoleMomentsCIS = requiresAllTransitionDipoleMomentsCIS; |
676 | 684 | } |
677 | 685 | |
686 | +vector<int>* Parameters::GetElectronicStateIndecesMullikenCIS() const{ | |
687 | + if(this->electronicStateIndecesMullikenCIS==NULL){ | |
688 | + stringstream ss; | |
689 | + ss << this->errorMessageGetElectronicStateIndecesMullikenCISNull; | |
690 | + throw MolDSException(ss.str()); | |
691 | + } | |
692 | + return this->electronicStateIndecesMullikenCIS; | |
693 | +} | |
694 | + | |
695 | +void Parameters::AddElectronicStateIndexMullikenCIS(int electronicStateIndex){ | |
696 | + if(this->electronicStateIndecesMullikenCIS==NULL){ | |
697 | + this->electronicStateIndecesMullikenCIS = new vector<int>; | |
698 | + } | |
699 | + this->electronicStateIndecesMullikenCIS->push_back(electronicStateIndex); | |
700 | +} | |
701 | + | |
702 | +bool Parameters::RequiresMullikenCIS() const{ | |
703 | + return (this->electronicStateIndecesMullikenCIS!=NULL && | |
704 | + 0<this->electronicStateIndecesMullikenCIS->size()); | |
705 | +} | |
706 | + | |
678 | 707 | // methods for Memory |
679 | 708 | double Parameters::GetLimitHeapMemory() const{ |
680 | 709 | return this->limitHeapMemory; |
@@ -115,28 +115,31 @@ public: | ||
115 | 115 | void SetRotatingEularAngles(double alpha, double beta, double gamma); |
116 | 116 | EularAngle GetRotatingEularAngles() const; |
117 | 117 | // CIS |
118 | - int GetActiveOccCIS() const; | |
119 | - void SetActiveOccCIS(int activeOccCIS); | |
120 | - int GetActiveVirCIS() const; | |
121 | - void SetActiveVirCIS(int activeOccCIS); | |
122 | - int GetNumberExcitedStatesCIS() const; | |
123 | - void SetNumberExcitedStatesCIS(int nStates); | |
124 | - bool RequiresCIS() const; | |
125 | - void SetRequiresCIS(bool requiresCIS); | |
126 | - bool IsDavidsonCIS() const; | |
127 | - void SetIsDavidsonCIS(bool isDavidsonCIS); | |
128 | - int GetMaxIterationsCIS() const; | |
129 | - void SetMaxIterationsCIS(int maxIterationsCIS); | |
130 | - int GetMaxDimensionsCIS() const; | |
131 | - void SetMaxDimensionsCIS(int maxDimensionsCIS); | |
132 | - double GetNormToleranceCIS() const; | |
133 | - void SetNormToleranceCIS(double normToleranceCIS); | |
134 | - int GetNumberPrintCoefficientsCIS() const; | |
135 | - void SetNumberPrintCoefficientsCIS(int numberPrintCoefficientsCIS); | |
136 | - bool RequiresExcitonEnergiesCIS() const; | |
137 | - void SetRequiresExcitonEnergiesCIS(bool requiresExcitonEnergiesCIS); | |
138 | - bool RequiresAllTransitionDipoleMomentsCIS() const; | |
139 | - void SetRequiresAllTransitionDipoleMomentsCIS(bool requiresAllTransitionDipoleMomentsCIS); | |
118 | + int GetActiveOccCIS() const; | |
119 | + void SetActiveOccCIS(int activeOccCIS); | |
120 | + int GetActiveVirCIS() const; | |
121 | + void SetActiveVirCIS(int activeOccCIS); | |
122 | + int GetNumberExcitedStatesCIS() const; | |
123 | + void SetNumberExcitedStatesCIS(int nStates); | |
124 | + bool RequiresCIS() const; | |
125 | + void SetRequiresCIS(bool requiresCIS); | |
126 | + bool IsDavidsonCIS() const; | |
127 | + void SetIsDavidsonCIS(bool isDavidsonCIS); | |
128 | + int GetMaxIterationsCIS() const; | |
129 | + void SetMaxIterationsCIS(int maxIterationsCIS); | |
130 | + int GetMaxDimensionsCIS() const; | |
131 | + void SetMaxDimensionsCIS(int maxDimensionsCIS); | |
132 | + double GetNormToleranceCIS() const; | |
133 | + void SetNormToleranceCIS(double normToleranceCIS); | |
134 | + int GetNumberPrintCoefficientsCIS() const; | |
135 | + void SetNumberPrintCoefficientsCIS(int numberPrintCoefficientsCIS); | |
136 | + bool RequiresExcitonEnergiesCIS() const; | |
137 | + void SetRequiresExcitonEnergiesCIS(bool requiresExcitonEnergiesCIS); | |
138 | + bool RequiresAllTransitionDipoleMomentsCIS() const; | |
139 | + void SetRequiresAllTransitionDipoleMomentsCIS(bool requiresAllTransitionDipoleMomentsCIS); | |
140 | + std::vector<int>* GetElectronicStateIndecesMullikenCIS() const; | |
141 | + void AddElectronicStateIndexMullikenCIS(int electronicStateIndex); | |
142 | + bool RequiresMullikenCIS() const; | |
140 | 143 | // Memory |
141 | 144 | double GetLimitHeapMemory() const; |
142 | 145 | void SetLimitHeapMemory(double limitHeap); |
@@ -214,6 +217,7 @@ private: | ||
214 | 217 | std::string errorMessageGetIndecesMOPlotNull; |
215 | 218 | std::string errorMessageGetIndecesHolePlotNull; |
216 | 219 | std::string errorMessageGetIndecesParticlePlotNull; |
220 | + std::string errorMessageGetElectronicStateIndecesMullikenCISNull; | |
217 | 221 | SimulationType currentSimulation; |
218 | 222 | TheoryType currentTheory; |
219 | 223 | // Physical constants |
@@ -268,17 +272,18 @@ private: | ||
268 | 272 | RotatingType rotatingType; |
269 | 273 | EularAngle rotatingEularAngles; |
270 | 274 | // CIS |
271 | - int activeOccCIS; | |
272 | - int activeVirCIS; | |
273 | - int numberExcitedStatesCIS; | |
274 | - int maxIterationsCIS; | |
275 | - int maxDimensionsCIS; | |
276 | - double normToleranceCIS; | |
277 | - bool requiresCIS; | |
278 | - bool isDavidsonCIS; | |
279 | - int numberPrintCoefficientsCIS; | |
280 | - bool requiresExcitonEnergiesCIS; | |
281 | - bool requiresAllTransitionDipoleMomentsCIS; | |
275 | + int activeOccCIS; | |
276 | + int activeVirCIS; | |
277 | + int numberExcitedStatesCIS; | |
278 | + int maxIterationsCIS; | |
279 | + int maxDimensionsCIS; | |
280 | + double normToleranceCIS; | |
281 | + bool requiresCIS; | |
282 | + bool isDavidsonCIS; | |
283 | + int numberPrintCoefficientsCIS; | |
284 | + bool requiresExcitonEnergiesCIS; | |
285 | + bool requiresAllTransitionDipoleMomentsCIS; | |
286 | + std::vector<int>* electronicStateIndecesMullikenCIS; | |
282 | 287 | // Memory |
283 | 288 | double limitHeapMemory; |
284 | 289 | // MD |
@@ -68,8 +68,10 @@ Cndo2::Cndo2(){ | ||
68 | 68 | this->matrixCISdimension = 0; |
69 | 69 | this->fockMatrix = NULL; |
70 | 70 | this->energiesMO = NULL; |
71 | - this->orbitalElectronPopulation = NULL; | |
72 | - this->atomicElectronPopulation = NULL; | |
71 | + this->orbitalElectronPopulation = NULL; | |
72 | + this->orbitalElectronPopulationCIS = NULL; | |
73 | + this->atomicElectronPopulation = NULL; | |
74 | + this->atomicElectronPopulationCIS = NULL; | |
73 | 75 | this->overlapAOs = NULL; |
74 | 76 | this->twoElecTwoCore = NULL; |
75 | 77 | this->cartesianMatrix = NULL; |
@@ -84,6 +84,8 @@ protected: | ||
84 | 84 | std::string messageStartSCF; |
85 | 85 | std::string messageDoneSCF; |
86 | 86 | std::string messageOmpElapsedTimeSCF; |
87 | + std::string messageMullikenAtoms; | |
88 | + std::string messageMullikenAtomsTitle; | |
87 | 89 | std::string messageUnitSec; |
88 | 90 | std::vector<MolDS_base::AtomType> enableAtomTypes; |
89 | 91 | MolDS_base::Molecule* molecule; |
@@ -94,7 +96,9 @@ protected: | ||
94 | 96 | double** fockMatrix; |
95 | 97 | double* energiesMO; |
96 | 98 | double** orbitalElectronPopulation; //P_{\mu\nu} of (2.50) in J. A. Pople book. |
99 | + double*** orbitalElectronPopulationCIS; | |
97 | 100 | double* atomicElectronPopulation; //P_{AB} of (3.21) in J. A. Pople book. |
101 | + double** atomicElectronPopulationCIS; | |
98 | 102 | double** overlapAOs; // overlap integral between AOs |
99 | 103 | double****** twoElecTwoCore; |
100 | 104 | double*** cartesianMatrix; // cartesian matrix represented by AOs |
@@ -250,8 +254,6 @@ private: | ||
250 | 254 | std::string messageDensityRMS; |
251 | 255 | std::string messageEnergyMO; |
252 | 256 | std::string messageEnergyMOTitle; |
253 | - std::string messageMullikenAtoms; | |
254 | - std::string messageMullikenAtomsTitle; | |
255 | 257 | std::string messageElecEnergy; |
256 | 258 | std::string messageNoteElecEnergy; |
257 | 259 | std::string messageNoteElecEnergyVdW; |
@@ -87,6 +87,15 @@ ZindoS::~ZindoS(){ | ||
87 | 87 | this->matrixForceElecStatesNum, |
88 | 88 | this->molecule->GetNumberAtoms(), |
89 | 89 | CartesianType_end); |
90 | + if(Parameters::GetInstance()->RequiresMullikenCIS()){ | |
91 | + MallocerFreer::GetInstance()->Free<double>(&this->orbitalElectronPopulationCIS, | |
92 | + Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS()->size(), | |
93 | + this->molecule->GetTotalNumberAOs(), | |
94 | + this->molecule->GetTotalNumberAOs()); | |
95 | + MallocerFreer::GetInstance()->Free<double>(&this->atomicElectronPopulationCIS, | |
96 | + Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS()->size(), | |
97 | + this->molecule->GetNumberAtoms()); | |
98 | + } | |
90 | 99 | //this->OutputLog("ZindoS deleted\n"); |
91 | 100 | } |
92 | 101 |
@@ -999,6 +1008,18 @@ void ZindoS::CalcCISProperties(){ | ||
999 | 1008 | this->energiesMO, |
1000 | 1009 | this->matrixCIS, |
1001 | 1010 | this->matrixCISdimension); |
1011 | + | |
1012 | + // orbital electron population | |
1013 | + this->CalcOrbitalElectronPopulationCIS(&this->orbitalElectronPopulationCIS, | |
1014 | + this->orbitalElectronPopulation, | |
1015 | + *this->molecule, | |
1016 | + this->fockMatrix, | |
1017 | + this->matrixCIS); | |
1018 | + | |
1019 | + // atomic electron population | |
1020 | + this->CalcAtomicElectronPopulationCIS(&this->atomicElectronPopulationCIS, | |
1021 | + this->orbitalElectronPopulationCIS, | |
1022 | + *this->molecule); | |
1002 | 1023 | } |
1003 | 1024 | |
1004 | 1025 | void ZindoS::CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments, |
@@ -1271,6 +1292,96 @@ void ZindoS::CalcFreeExcitonEnergies(double** freeExcitonEnergiesCIS, | ||
1271 | 1292 | } |
1272 | 1293 | } |
1273 | 1294 | |
1295 | +void ZindoS::CalcOrbitalElectronPopulationCIS(double**** orbitalElectronPopulationCIS, | |
1296 | + double const* const* orbitalElectronPopulation, | |
1297 | + const MolDS_base::Molecule& molecule, | |
1298 | + double const* const* fockMatrix, | |
1299 | + double const* const* matrixCIS) const{ | |
1300 | + if(!Parameters::GetInstance()->RequiresMullikenCIS()){ | |
1301 | + return; | |
1302 | + } | |
1303 | + vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS(); | |
1304 | + // malloc or initialize free exciton energies | |
1305 | + if(*orbitalElectronPopulationCIS == NULL){ | |
1306 | + MallocerFreer::GetInstance()->Malloc<double>(orbitalElectronPopulationCIS, | |
1307 | + elecStates->size(), | |
1308 | + molecule.GetTotalNumberAOs(), | |
1309 | + molecule.GetTotalNumberAOs()); | |
1310 | + } | |
1311 | + else{ | |
1312 | + MallocerFreer::GetInstance()->Initialize<double>(*orbitalElectronPopulationCIS, | |
1313 | + elecStates->size(), | |
1314 | + molecule.GetTotalNumberAOs(), | |
1315 | + molecule.GetTotalNumberAOs()); | |
1316 | + } | |
1317 | + // clac orbital electron population | |
1318 | + int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2; | |
1319 | + int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS(); | |
1320 | + int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS(); | |
1321 | + for(int k=0; k<elecStates->size(); k++){ | |
1322 | + int excitedStateIndex = (*elecStates)[k]-1; | |
1323 | + for(int mu=0; mu<molecule.GetTotalNumberAOs(); mu++){ | |
1324 | + for(int nu=0; nu<molecule.GetTotalNumberAOs(); nu++){ | |
1325 | + double value = orbitalElectronPopulation[mu][nu]; | |
1326 | + for(int moI=0; moI<numberActiveOcc; moI++){ | |
1327 | + for(int moA=numberOcc; moA<numberOcc+numberActiveVir; moA++){ | |
1328 | + int slaterDeterminantIndex = this->GetSlaterDeterminantIndex(moI,moA); | |
1329 | + value += pow(matrixCIS[excitedStateIndex][slaterDeterminantIndex],2.0) | |
1330 | + *(-fockMatrix[moI][mu]*fockMatrix[moI][nu] | |
1331 | + +fockMatrix[moA][mu]*fockMatrix[moA][nu]); | |
1332 | + double tmpVal1=0.0; | |
1333 | + for(int moB=numberOcc; moB<numberOcc+numberActiveVir; moB++){ | |
1334 | + if(moB==moA) continue; | |
1335 | + int tmpSDIndex = this->GetSlaterDeterminantIndex(moI,moB); | |
1336 | + tmpVal1 += matrixCIS[excitedStateIndex][tmpSDIndex]*fockMatrix[moB][nu]; | |
1337 | + } | |
1338 | + double tmpVal2=0.0; | |
1339 | + for(int moJ=0; moJ<numberActiveOcc; moJ++){ | |
1340 | + if(moJ==moI) continue; | |
1341 | + int tmpSDIndex = this->GetSlaterDeterminantIndex(moJ,moA); | |
1342 | + tmpVal2 += matrixCIS[excitedStateIndex][tmpSDIndex]*fockMatrix[moJ][mu]; | |
1343 | + } | |
1344 | + value += matrixCIS[excitedStateIndex][slaterDeterminantIndex] | |
1345 | + *(fockMatrix[moA][mu]*tmpVal1 + fockMatrix[moI][nu]*tmpVal2); | |
1346 | + } | |
1347 | + } | |
1348 | + (*orbitalElectronPopulationCIS)[k][mu][nu] = value; | |
1349 | + } | |
1350 | + } | |
1351 | + } | |
1352 | +} | |
1353 | + | |
1354 | +void ZindoS::CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationCIS, | |
1355 | + double const* const* const* orbitalElectronPopulationCIS, | |
1356 | + const Molecule& molecule) const{ | |
1357 | + if(!Parameters::GetInstance()->RequiresMullikenCIS()){ | |
1358 | + return; | |
1359 | + } | |
1360 | + int totalNumberAtoms = molecule.GetNumberAtoms(); | |
1361 | + vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS(); | |
1362 | + // malloc or initialize free exciton energies | |
1363 | + if(*atomicElectronPopulationCIS == NULL){ | |
1364 | + MallocerFreer::GetInstance()->Malloc<double>(atomicElectronPopulationCIS, | |
1365 | + elecStates->size(), | |
1366 | + totalNumberAtoms); | |
1367 | + } | |
1368 | + else{ | |
1369 | + MallocerFreer::GetInstance()->Initialize<double>(*atomicElectronPopulationCIS, | |
1370 | + elecStates->size(), | |
1371 | + totalNumberAtoms); | |
1372 | + } | |
1373 | + // clac atomic electron population | |
1374 | + for(int k=0; k<elecStates->size(); k++){ | |
1375 | + for(int a=0; a<totalNumberAtoms; a++){ | |
1376 | + int firstAOIndex = molecule.GetAtom(a)->GetFirstAOIndex(); | |
1377 | + int numberAOs = molecule.GetAtom(a)->GetValenceSize(); | |
1378 | + for(int i=firstAOIndex; i<firstAOIndex+numberAOs; i++){ | |
1379 | + (*atomicElectronPopulationCIS)[k][a] += orbitalElectronPopulationCIS[k][i][i]; | |
1380 | + } | |
1381 | + } | |
1382 | + } | |
1383 | +} | |
1384 | + | |
1274 | 1385 | void ZindoS::OutputCISResults() const{ |
1275 | 1386 | int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS(); |
1276 | 1387 | int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS(); |
@@ -1303,6 +1414,9 @@ void ZindoS::OutputCISResults() const{ | ||
1303 | 1414 | // output transition dipole moment |
1304 | 1415 | this->OutputCISTransitionDipole(); |
1305 | 1416 | |
1417 | + // output mulliken population | |
1418 | + this->OutputCISMulliken(); | |
1419 | + | |
1306 | 1420 | // output exciton energies |
1307 | 1421 | if(Parameters::GetInstance()->RequiresExcitonEnergiesCIS()){ |
1308 | 1422 | this->OutputLog(this->messageExcitonEnergiesCIS); |
@@ -1427,8 +1541,29 @@ void ZindoS::OutputCISTransitionDipole() const{ | ||
1427 | 1541 | } |
1428 | 1542 | } |
1429 | 1543 | this->OutputLog("\n"); |
1544 | +} | |
1430 | 1545 | |
1546 | +void ZindoS::OutputCISMulliken() const{ | |
1547 | + if(!Parameters::GetInstance()->RequiresMullikenCIS()){ | |
1548 | + return; | |
1549 | + } | |
1550 | + int totalNumberAtoms = this->molecule->GetNumberAtoms(); | |
1551 | + this->OutputLog(this->messageMullikenAtomsTitle); | |
1552 | + vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS(); | |
1553 | + for(int k=0; k<elecStates->size(); k++){ | |
1554 | + for(int a=0; a<totalNumberAtoms; a++){ | |
1555 | + const Atom& atom = *this->molecule->GetAtom(a); | |
1556 | + this->OutputLog(boost::format("%s\t%d\t%d\t%s\t%e\t%e\n") % this->messageMullikenAtoms | |
1557 | + % (*elecStates)[k] | |
1558 | + % a | |
1559 | + % AtomTypeStr(atom.GetAtomType()) | |
1560 | + % atom.GetCoreCharge() | |
1561 | + % (atom.GetCoreCharge()-this->atomicElectronPopulationCIS[k][a])); | |
1562 | + } | |
1563 | + this->OutputLog("\n"); | |
1564 | + } | |
1431 | 1565 | } |
1566 | + | |
1432 | 1567 | void ZindoS::SortCISEigenVectorCoefficients(vector<CISEigenVectorCoefficient>* cisEigenVectorCoefficients, |
1433 | 1568 | double* cisEigenVector) const{ |
1434 | 1569 | for(int l=0; l<this->matrixCISdimension; l++){ |
@@ -138,11 +138,20 @@ private: | ||
138 | 138 | void DoCISDavidson(); |
139 | 139 | void OutputCISDipole() const; |
140 | 140 | void OutputCISTransitionDipole() const; |
141 | + void OutputCISMulliken() const; | |
141 | 142 | void CalcFreeExcitonEnergies(double** freeExcitonEnergiesCIS, |
142 | 143 | const MolDS_base::Molecule& molecule, |
143 | 144 | double const* energiesMO, |
144 | 145 | double const* const* matrixCIS, |
145 | 146 | int matrixCISdimension) const; |
147 | + void CalcOrbitalElectronPopulationCIS(double**** orbitalElectronPopulationCIS, | |
148 | + double const* const* orbitalElectronPopulation, | |
149 | + const MolDS_base::Molecule& molecule, | |
150 | + double const* const* fockMatrix, | |
151 | + double const* const* matrixCIS) const; | |
152 | + void CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationCIS, | |
153 | + double const* const* const* orbitalElectronPopulationCIS, | |
154 | + const MolDS_base::Molecule& molecule) const; | |
146 | 155 | void CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments, |
147 | 156 | double const* const* fockMatrix, |
148 | 157 | double const* const* matrixCIS, |