votca 2024-dev
Loading...
Searching...
No Matches
bsecoupling.cc
Go to the documentation of this file.
1/*
2 * Copyright 2009-2020 The VOTCA Development Team
3 * (http://www.votca.org)
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License")
6 *
7 * You may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 *
18 */
19
20// Third party includes
21#include <boost/format.hpp>
22
23// VOTCA includes
24#include <votca/tools/constants.h>
25
26// Local VOTCA includes
27#include "votca/xtp/aomatrix.h"
28#include "votca/xtp/bse.h"
29#include "votca/xtp/bse_operator.h"
30#include "votca/xtp/bsecoupling.h"
31
32namespace votca {
33namespace xtp {
34using namespace std;
35using boost::format;
36using namespace tools;
37
38void BSECoupling::Initialize(Property& options) {
39 string spintype = options.get("spin").as<std::string>();
40 if (spintype == "all") {
41 doSinglets_ = true;
42 doTriplets_ = true;
43 } else if (spintype == "triplet") {
44 doTriplets_ = true;
45 doSinglets_ = false;
46 } else if (spintype == "singlet") {
47 doSinglets_ = true;
48 doTriplets_ = false;
49 } else {
50 throw std::runtime_error(
51 (boost::format(
52 "Choice % for type not known. Available singlet,triplet,all") %
53 spintype)
54 .str());
55 }
56 output_perturbation_ = options.get("use_perturbation").as<bool>();
57
58 levA_ = options.get("moleculeA.states").as<Index>();
59 levB_ = options.get("moleculeB.states").as<Index>();
60 occA_ = options.get("moleculeA.occLevels").as<Index>();
61 occB_ = options.get("moleculeB.occLevels").as<Index>();
62 unoccA_ = options.get("moleculeA.unoccLevels").as<Index>();
63 unoccB_ = options.get("moleculeB.unoccLevels").as<Index>();
64}
65
66void BSECoupling::WriteToProperty(Property& summary, const QMState& stateA,
67 const QMState& stateB) const {
68 Property& coupling_summary = summary.add("coupling", "");
69 double JAB_pert = 0;
70 double JAB_diag = 0;
71 if (stateA.Type() == QMStateType::Singlet) {
72 JAB_pert =
73 getSingletCouplingElement(stateA.StateIdx(), stateB.StateIdx(), 0);
74 JAB_diag =
75 getSingletCouplingElement(stateA.StateIdx(), stateB.StateIdx(), 1);
76 } else if (stateA.Type() == QMStateType::Triplet) {
77 JAB_pert =
78 getTripletCouplingElement(stateA.StateIdx(), stateB.StateIdx(), 0);
79 JAB_diag =
80 getTripletCouplingElement(stateA.StateIdx(), stateB.StateIdx(), 1);
81 }
82 coupling_summary.setAttribute("stateA", stateA.ToString());
83 coupling_summary.setAttribute("stateB", stateB.ToString());
84 coupling_summary.setAttribute("j_pert", (format("%1$1.6e") % JAB_pert).str());
85 coupling_summary.setAttribute("j_diag", (format("%1$1.6e") % JAB_diag).str());
86}
87
88void BSECoupling::Addoutput(Property& type_summary, const Orbitals&,
89 const Orbitals&) const {
90 tools::Property& bsecoupling = type_summary.add(Identify(), "");
91 string algorithm = "j_diag";
92 if (output_perturbation_) {
93 algorithm = "j_pert";
94 }
95 if (doSinglets_) {
96 QMStateType singlet = QMStateType(QMStateType::Singlet);
97 Property& singlet_summary = bsecoupling.add(singlet.ToLongString(), "");
98 singlet_summary.setAttribute("algorithm", algorithm);
99 for (Index stateA = 0; stateA < levA_; ++stateA) {
100 QMState qmstateA = QMState(singlet, stateA, false);
101 for (Index stateB = 0; stateB < levB_; ++stateB) {
102 QMState qmstateB = QMState(singlet, stateB, false);
103 WriteToProperty(singlet_summary, qmstateA, qmstateB);
104 }
105 }
106 }
107
108 if (doTriplets_) {
109 QMStateType triplet = QMStateType(QMStateType::Triplet);
110 Property& triplet_summary = bsecoupling.add(triplet.ToLongString(), "");
111 triplet_summary.setAttribute("algorithm", algorithm);
112 for (Index stateA = 0; stateA < levA_; ++stateA) {
113 QMState qmstateA = QMState(triplet, stateA, false);
114 for (Index stateB = 0; stateB < levB_; ++stateB) {
115 QMState qmstateB = QMState(triplet, stateB, false);
116 WriteToProperty(triplet_summary, qmstateA, qmstateB);
117 }
118 }
119 }
120}
121
122double BSECoupling::getSingletCouplingElement(Index levelA, Index levelB,
123 Index methodindex) const {
124 return JAB_singlet[methodindex](levelA, levelB + levA_) *
125 votca::tools::conv::hrt2ev;
126}
127
128double BSECoupling::getTripletCouplingElement(Index levelA, Index levelB,
129 Index methodindex) const {
130 return JAB_triplet[methodindex](levelA, levelB + levA_) *
131 votca::tools::conv::hrt2ev;
132}
133
134Eigen::MatrixXd BSECoupling::SetupCTStates(Index bseA_vtotal, Index bseB_vtotal,
135 Index bseAB_vtotal,
136 Index bseAB_ctotal,
137 const Eigen::MatrixXd& A_AB,
138 const Eigen::MatrixXd& B_AB) const {
139
140 Index noAB = occA_ * unoccB_;
141 Index noBA = unoccA_ * occB_;
142 Index bseAB_size = bseAB_vtotal * bseAB_ctotal;
143 Eigen::MatrixXd CTstates = Eigen::MatrixXd::Zero(bseAB_size, noAB + noBA);
144
145 auto A_occ = A_AB.middleCols(bseA_vtotal - occA_, occA_);
146 auto A_unocc = A_AB.middleCols(bseA_vtotal, unoccA_);
147 auto B_occ = B_AB.middleCols(bseB_vtotal - occB_, occB_);
148 auto B_unocc = B_AB.middleCols(bseB_vtotal, unoccB_);
149
150 const Eigen::MatrixXd A_occ_occ = A_occ.topRows(bseAB_vtotal);
151 const Eigen::MatrixXd B_unocc_unocc = B_unocc.bottomRows(bseAB_ctotal);
152
153 // notation AB is CT states with A+B-, BA is the counterpart
154 // Setting up CT-states:
155 XTP_LOG(Log::error, *pLog_)
156 << TimeStamp() << " Setting up CT-states" << flush;
157
158 // Number of A+B- states
159#pragma omp parallel for
160 for (Index a_occ = 0; a_occ < occA_; a_occ++) {
161 for (Index b_unocc = 0; b_unocc < unoccB_; b_unocc++) {
162 Index index = a_occ * unoccB_ + b_unocc;
163 Eigen::MatrixXd Coeff =
164 B_unocc_unocc.col(b_unocc) * A_occ_occ.col(a_occ).transpose();
165 CTstates.col(index) =
166 Eigen::Map<Eigen::VectorXd>(Coeff.data(), bseAB_size);
167 }
168 }
169 XTP_LOG(Log::error, *pLog_)
170 << TimeStamp() << " " << noBA << " CT states A+B- created" << flush;
171
172 const Eigen::MatrixXd A_unocc_unocc = A_unocc.bottomRows(bseAB_ctotal);
173 const Eigen::MatrixXd B_occ_occ = B_occ.topRows(bseAB_vtotal);
174
175#pragma omp parallel for
176 for (Index b_occ = 0; b_occ < occB_; b_occ++) {
177 for (Index a_unocc = 0; a_unocc < unoccA_; a_unocc++) {
178 Index index = b_occ * unoccA_ + a_unocc + noAB;
179 Eigen::MatrixXd Coeff =
180 A_unocc_unocc.col(a_unocc) * B_occ_occ.col(b_occ).transpose();
181 CTstates.col(index) =
182 Eigen::Map<Eigen::VectorXd>(Coeff.data(), bseAB_size);
183 }
184 }
185 XTP_LOG(Log::error, *pLog_)
186 << TimeStamp() << " " << noBA << " CT states A-B+ created" << flush;
187 return CTstates;
188}
189
190Eigen::MatrixXd BSECoupling::ProjectFrenkelExcitons(
191 const Eigen::MatrixXd& BSE_Coeffs, const Eigen::MatrixXd& X_AB,
192 Index bseX_vtotal, Index bseX_ctotal, Index bseAB_vtotal,
193 Index bseAB_ctotal) const {
194 Index bseAB_size = bseAB_vtotal * bseAB_ctotal;
195 auto X_occ = X_AB.leftCols(bseX_vtotal);
196 auto X_unocc = X_AB.rightCols(bseX_ctotal);
197 const Eigen::MatrixXd X_occ_occ = X_occ.topRows(bseAB_vtotal);
198 const Eigen::MatrixXd X_unocc_unocc = X_unocc.bottomRows(bseAB_ctotal);
199 Eigen::MatrixXd result = Eigen::MatrixXd::Zero(bseAB_size, BSE_Coeffs.cols());
200 // no pragma here because often we will only have one Coeff
201 for (Index i = 0; i < BSE_Coeffs.cols(); i++) {
202 Eigen::VectorXd coeff = BSE_Coeffs.col(i);
203 Eigen::Map<Eigen::MatrixXd> coeffmatrix =
204 Eigen::Map<Eigen::MatrixXd>(coeff.data(), bseX_ctotal, bseX_vtotal);
205 Eigen::MatrixXd proj = X_unocc_unocc * coeffmatrix * X_occ_occ.transpose();
206 result.col(i) = Eigen::Map<Eigen::VectorXd>(proj.data(), proj.size());
207 }
208 return result;
209}
210
211int GetSign(double value) {
212 if (value < 0) {
213 return -1;
214 } else if (value > 0) {
215 return 1;
216 }
217 return 0;
218}
219
220void BSECoupling::CalculateCouplings(const Orbitals& orbitalsA,
221 const Orbitals& orbitalsB,
222 const Orbitals& orbitalsAB) {
223 XTP_LOG(Log::error, *pLog_)
224 << TimeStamp() << " Calculating exciton couplings" << flush;
225 // set the parallelization
226 XTP_LOG(Log::error, *pLog_) << TimeStamp() << " Using "
227 << OPENMP::getMaxThreads() << " threads" << flush;
228
229 CheckAtomCoordinates(orbitalsA, orbitalsB, orbitalsAB);
230
231 // constructing the direct product orbA x orbB
232 Index basisB = orbitalsB.getBasisSetSize();
233 Index basisA = orbitalsA.getBasisSetSize();
234 if ((basisA == 0) || (basisB == 0)) {
235 throw std::runtime_error("Basis set size is not stored in monomers");
236 }
237 // get exciton information of molecule A
238 Index bseA_cmax = orbitalsA.getBSEcmax();
239 Index bseA_cmin = orbitalsA.getBSEcmin();
240 Index bseA_vmax = orbitalsA.getBSEvmax();
241 Index bseA_vmin = orbitalsA.getBSEvmin();
242 Index bseA_vtotal = bseA_vmax - bseA_vmin + 1;
243 Index bseA_ctotal = bseA_cmax - bseA_cmin + 1;
244 Index bseA_total = bseA_vtotal + bseA_ctotal;
245 Index bseA_size = bseA_vtotal * bseA_ctotal;
246 Index bseA_singlet_exc = orbitalsA.BSESinglets().eigenvectors().cols();
247 Index bseA_triplet_exc = orbitalsA.BSETriplets().eigenvectors().cols();
248
249 XTP_LOG(Log::error, *pLog_)
250 << TimeStamp() << " molecule A has " << bseA_singlet_exc
251 << " singlet excitons with dimension " << bseA_size << flush;
252 XTP_LOG(Log::error, *pLog_)
253 << TimeStamp() << " molecule A has " << bseA_triplet_exc
254 << " triplet excitons with dimension " << bseA_size << flush;
255
256 // get exciton information of molecule B
257 Index bseB_cmax = orbitalsB.getBSEcmax();
258 Index bseB_cmin = orbitalsB.getBSEcmin();
259 Index bseB_vmax = orbitalsB.getBSEvmax();
260 Index bseB_vmin = orbitalsB.getBSEvmin();
261 Index bseB_vtotal = bseB_vmax - bseB_vmin + 1;
262 Index bseB_ctotal = bseB_cmax - bseB_cmin + 1;
263 Index bseB_total = bseB_vtotal + bseB_ctotal;
264 Index bseB_size = bseB_vtotal * bseB_ctotal;
265 Index bseB_singlet_exc = orbitalsB.BSESinglets().eigenvectors().cols();
266 Index bseB_triplet_exc = orbitalsB.BSETriplets().eigenvectors().cols();
267
268 XTP_LOG(Log::error, *pLog_)
269 << TimeStamp() << " molecule B has " << bseB_singlet_exc
270 << " singlet excitons with dimension " << bseB_size << flush;
271 XTP_LOG(Log::error, *pLog_)
272 << TimeStamp() << " molecule B has " << bseB_triplet_exc
273 << " triplet excitons with dimension " << bseB_size << flush;
274
275 if (doSinglets_) {
276 if (levA_ > bseA_singlet_exc) {
277 XTP_LOG(Log::error, *pLog_) << TimeStamp()
278 << " Number of excitons you want is greater "
279 "than stored for molecule "
280 "A. Setting to max number available"
281 << flush;
282 levA_ = bseA_singlet_exc;
283 }
284 if (levB_ > bseB_singlet_exc) {
285 XTP_LOG(Log::error, *pLog_) << TimeStamp()
286 << " Number of excitons you want is greater "
287 "than stored for molecule "
288 "B. Setting to max number available"
289 << flush;
290 levB_ = bseB_singlet_exc;
291 }
292 }
293 if (doTriplets_) {
294 if (levA_ > bseA_triplet_exc) {
295 XTP_LOG(Log::error, *pLog_)
296 << TimeStamp()
297 << " Number of Frenkel states you want is greater than stored for "
298 "molecule A. Setting to max number available"
299 << flush;
300 levA_ = bseA_triplet_exc;
301 }
302 if (levB_ > bseB_triplet_exc) {
303 XTP_LOG(Log::error, *pLog_)
304 << TimeStamp()
305 << " Number of Frenkel states you want is greater than stored for "
306 "molecule B. Setting to max number available"
307 << flush;
308 levB_ = bseB_triplet_exc;
309 }
310 }
311 if (unoccA_ > bseA_ctotal || unoccA_ < 0) {
312 XTP_LOG(Log::error, *pLog_)
313 << TimeStamp()
314 << " Number of unoccupied orbitals in molecule A for CT creation "
315 "exceeds number of KS-orbitals in BSE"
316 << flush;
317 unoccA_ = bseA_ctotal;
318 }
319 if (unoccB_ > bseB_ctotal || unoccB_ < 0) {
320 XTP_LOG(Log::error, *pLog_)
321 << TimeStamp()
322 << " Number of unoccupied orbitals in molecule B for CT creation "
323 "exceeds number of KS-orbitals in BSE"
324 << flush;
325 unoccB_ = bseB_ctotal;
326 }
327
328 if (occA_ > bseA_vtotal || occA_ < 0) {
329 XTP_LOG(Log::error, *pLog_)
330 << TimeStamp()
331 << " Number of occupied orbitals in molecule A for CT creation "
332 "exceeds number of KS-orbitals in BSE"
333 << flush;
334 occA_ = bseA_vtotal;
335 }
336 if (occB_ > bseB_vtotal || occB_ < 0) {
337 XTP_LOG(Log::error, *pLog_)
338 << TimeStamp()
339 << " Number of occupied orbitals in molecule B for CT creation "
340 "exceeds number of KS-orbitals in BSE"
341 << flush;
342 occB_ = bseB_vtotal;
343 }
344
345 // get exciton information of pair AB
346 Index bseAB_cmax = orbitalsAB.getBSEcmax();
347 Index bseAB_cmin = orbitalsAB.getBSEcmin();
348 Index bseAB_vmax = orbitalsAB.getBSEvmax();
349 Index bseAB_vmin = orbitalsAB.getBSEvmin();
350 Index bseAB_vtotal = bseAB_vmax - bseAB_vmin + 1;
351 Index bseAB_ctotal = bseAB_cmax - bseAB_cmin + 1;
352 Index bseAB_total = bseAB_vtotal + bseAB_ctotal;
353 Index bseAB_size = bseAB_vtotal * bseAB_ctotal;
354
355 // DFT levels of monomers can be reduced to those used in BSE
356 XTP_LOG(Log::error, *pLog_)
357 << TimeStamp() << " levels used for BSE of molA: " << bseA_vmin
358 << " to " << bseA_cmax << " total: " << bseA_total << flush;
359 XTP_LOG(Log::error, *pLog_)
360 << TimeStamp() << " levels used for BSE of molB: " << bseB_vmin
361 << " to " << bseB_cmax << " total: " << bseB_total << flush;
362 XTP_LOG(Log::error, *pLog_)
363 << TimeStamp() << " levels used for BSE of dimer AB: " << bseAB_vmin
364 << " to " << bseAB_cmax << " total: " << bseAB_total << flush;
365
366 Eigen::MatrixXd MOsA =
367 orbitalsA.MOs().eigenvectors().middleCols(bseA_vmin, bseA_total);
368 Eigen::MatrixXd MOsB =
369 orbitalsB.MOs().eigenvectors().middleCols(bseB_vmin, bseB_total);
370 Eigen::MatrixXd MOsAB =
371 orbitalsAB.MOs().eigenvectors().middleCols(bseAB_vmin, bseAB_total);
372
373 XTP_LOG(Log::info, *pLog_)
374 << TimeStamp() << " Calculating overlap matrix for basisset: "
375 << orbitalsAB.getDFTbasisName() << flush;
376
377 Eigen::MatrixXd overlap = CalculateOverlapMatrix(orbitalsAB) * MOsAB;
378 XTP_LOG(Log::info, *pLog_)
379 << TimeStamp() << " Projecting monomers onto dimer orbitals" << flush;
380
381 Eigen::MatrixXd A_AB = overlap.topRows(basisA).transpose() * MOsA;
382 Eigen::MatrixXd B_AB = overlap.bottomRows(basisB).transpose() * MOsB;
383 Eigen::VectorXd mag_A = A_AB.colwise().squaredNorm();
384 if (mag_A.any() < 0.95) {
385 XTP_LOG(Log::error, *pLog_)
386 << "\nWarning: "
387 << "Projection of orbitals of monomer A on dimer is insufficient,mag="
388 << mag_A.minCoeff() << flush;
389 }
390 Eigen::VectorXd mag_B = B_AB.colwise().squaredNorm();
391 if (mag_B.any() < 0.95) {
392 XTP_LOG(Log::error, *pLog_)
393 << "\nWarning: "
394 << "Projection of orbitals of monomer B on dimer is insufficient,mag="
395 << mag_B.minCoeff() << flush;
396 }
397
398 AOBasis dftbasis = orbitalsAB.getDftBasis();
399 AOBasis auxbasis = orbitalsAB.getAuxBasis();
400 TCMatrix_gwbse Mmn;
401 // rpamin here, because RPA needs till rpamin
402 Mmn.Initialize(auxbasis.AOBasisSize(), orbitalsAB.getRPAmin(),
403 orbitalsAB.getGWAmax(), orbitalsAB.getRPAmin(),
404 orbitalsAB.getRPAmax());
405 Mmn.Fill(auxbasis, dftbasis, orbitalsAB.MOs().eigenvectors());
406
407 const Eigen::MatrixXd& qpcoeff = orbitalsAB.QPdiag().eigenvectors();
408 Eigen::MatrixXd Hqp = qpcoeff *
409 orbitalsAB.QPdiag().eigenvalues().asDiagonal() *
410 qpcoeff.transpose();
411
412 BSE::options opt;
413 opt.cmax = orbitalsAB.getBSEcmax();
414 opt.homo = orbitalsAB.getHomo();
415 opt.qpmin = orbitalsAB.getGWAmin();
416 opt.qpmax = orbitalsAB.getGWAmax();
417 opt.rpamax = orbitalsAB.getRPAmax();
418 opt.rpamin = orbitalsAB.getRPAmin();
419 opt.useTDA = true;
420 opt.vmin = orbitalsAB.getBSEvmin();
421 opt.use_Hqp_offdiag = orbitalsAB.GetFlagUseHqpOffdiag();
422 BSE bse(*pLog_, Mmn);
423 bse.configure(opt, orbitalsAB.RPAInputEnergies(), Hqp);
424 XTP_LOG(Log::error, *pLog_) << TimeStamp() << " Setup BSE operator" << flush;
425
426 // now the different spin types
427 if (doSinglets_) {
428 XTP_LOG(Log::error, *pLog_)
429 << TimeStamp() << " Evaluating singlets" << flush;
430 XTP_LOG(Log::error, *pLog_)
431 << TimeStamp() << " Setup Hamiltonian" << flush;
432 Eigen::MatrixXd FE_AB = Eigen::MatrixXd::Zero(bseAB_size, levA_ + levB_);
433 const Eigen::MatrixXd bseA =
434 orbitalsA.BSESinglets().eigenvectors().leftCols(levA_);
435 FE_AB.leftCols(levA_) = ProjectFrenkelExcitons(
436 bseA, A_AB, bseA_vtotal, bseA_ctotal, bseAB_vtotal, bseAB_ctotal);
437 const Eigen::MatrixXd bseB =
438 orbitalsB.BSESinglets().eigenvectors().leftCols(levB_);
439 FE_AB.rightCols(levB_) = ProjectFrenkelExcitons(
440 bseB, B_AB, bseB_vtotal, bseB_ctotal, bseAB_vtotal, bseAB_ctotal);
441 Eigen::MatrixXd CTStates = SetupCTStates(
442 bseA_vtotal, bseB_vtotal, bseAB_vtotal, bseAB_ctotal, A_AB, B_AB);
443 JAB_singlet =
444 ProjectExcitons(FE_AB, CTStates, bse.getSingletOperator_TDA());
445 XTP_LOG(Log::error, *pLog_)
446 << TimeStamp() << " calculated singlet couplings " << flush;
447 }
448 if (doTriplets_) {
449 XTP_LOG(Log::error, *pLog_)
450 << TimeStamp() << " Evaluating triplets" << flush;
451 Eigen::MatrixXd FE_AB = Eigen::MatrixXd::Zero(bseAB_size, levA_ + levB_);
452 const Eigen::MatrixXd bseA =
453 orbitalsA.BSETriplets().eigenvectors().leftCols(levA_);
454 FE_AB.leftCols(levA_) = ProjectFrenkelExcitons(
455 bseA, A_AB, bseA_vtotal, bseA_ctotal, bseAB_vtotal, bseAB_ctotal);
456 const Eigen::MatrixXd bseB =
457 orbitalsB.BSETriplets().eigenvectors().leftCols(levB_);
458 FE_AB.rightCols(levB_) = ProjectFrenkelExcitons(
459 bseB, B_AB, bseB_vtotal, bseB_ctotal, bseAB_vtotal, bseAB_ctotal);
460 Eigen::MatrixXd CTStates = SetupCTStates(
461 bseA_vtotal, bseB_vtotal, bseAB_vtotal, bseAB_ctotal, A_AB, B_AB);
462 JAB_triplet =
463 ProjectExcitons(FE_AB, CTStates, bse.getTripletOperator_TDA());
464 XTP_LOG(Log::error, *pLog_)
465 << TimeStamp() << " calculated triplet couplings " << flush;
466 }
467 XTP_LOG(Log::error, *pLog_)
468 << TimeStamp() << " Done with exciton couplings" << flush;
469}
470
471Eigen::MatrixXd BSECoupling::OrthogonalizeCTs(Eigen::MatrixXd& FE_AB,
472 Eigen::MatrixXd& CTStates) const {
473 Index ct = CTStates.cols();
474
475 if (ct > 0) {
476 XTP_LOG(Log::error, *pLog_)
477 << TimeStamp() << " Orthogonalizing CT-states with respect to FE-states"
478 << flush;
479 Eigen::MatrixXd correction = FE_AB * (FE_AB.transpose() * CTStates);
480 CTStates -= correction;
481
482 // normalize
483 Eigen::VectorXd norm = CTStates.colwise().norm();
484 for (Index i = 0; i < CTStates.cols(); i++) {
485 CTStates.col(i) /= norm(i);
486 }
487 Index minstateindex = 0;
488 double minnorm = norm.minCoeff(&minstateindex);
489 if (minnorm < 0.95) {
490 XTP_LOG(Log::error, *pLog_)
491 << TimeStamp() << " WARNING: CT-state " << minstateindex
492 << " norm is only " << minnorm << flush;
493 }
494 }
495 Index bse_exc = levA_ + levB_;
496
497 Index bseAB_size = CTStates.rows();
498 Eigen::MatrixXd projection(bseAB_size, bse_exc + ct);
499 XTP_LOG(Log::info, *pLog_)
500 << TimeStamp() << " merging projections into one vector " << flush;
501 projection.leftCols(bse_exc) = FE_AB;
502 FE_AB.resize(0, 0);
503 if (ct > 0) {
504 projection.rightCols(ct) = CTStates;
505 }
506 CTStates.resize(0, 0);
507 return projection;
508}
509
510template <class BSE_OPERATOR>
511Eigen::MatrixXd BSECoupling::CalcJ_dimer(BSE_OPERATOR& H,
512 Eigen::MatrixXd& projection) const {
513
514 XTP_LOG(Log::info, *pLog_)
515 << TimeStamp() << " Setting up coupling matrix size "
516 << projection.cols() << flush;
517 // matrix J_
518 // E_A J_AB J_A_ABCT J_A_BACT
519 // J_BA E_B J_B_ABCT J_B_BACT
520 // J_ABCT_A J_ABCT_B E_ABCT J_ABCT_BACT
521 // J_BACT_A J_BACT_B J_BACT_ABCT E_BACT
522
523 // this only works for hermitian/symmetric H so only in TDA
524
525 Eigen::MatrixXd J_dimer = projection.transpose() * (H * projection).eval();
526
527 XTP_LOG(Log::info, *pLog_)
528 << TimeStamp() << " Setting up overlap matrix size "
529 << projection.cols() << flush;
530 Eigen::MatrixXd S_dimer = projection.transpose() * projection;
531
532 projection.resize(0, 0);
533 if (projection.cols()) {
534 XTP_LOG(Log::debug, *pLog_)
535 << "---------------------------------------" << flush;
536 XTP_LOG(Log::debug, *pLog_) << " J_dimer_[Ryd]" << flush;
537
538 XTP_LOG(Log::debug, *pLog_) << J_dimer << flush;
539 XTP_LOG(Log::debug, *pLog_) << " S_dimer_" << flush;
540
541 XTP_LOG(Log::debug, *pLog_) << S_dimer << flush;
542 XTP_LOG(Log::debug, *pLog_)
543 << "---------------------------------------" << flush;
544 }
545
546 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(S_dimer);
547 Eigen::MatrixXd Sm1 = es.operatorInverseSqrt();
548 Eigen::MatrixXd J_ortho = Sm1 * J_dimer * Sm1;
549
550 if (projection.cols()) {
551 XTP_LOG(Log::debug, *pLog_)
552 << "---------------------------------------" << flush;
553 XTP_LOG(Log::debug, *pLog_) << " J_ortho_[Ryd]" << flush;
554 XTP_LOG(Log::debug, *pLog_) << J_ortho << flush;
555 XTP_LOG(Log::debug, *pLog_) << " S_-1/2" << flush;
556 XTP_LOG(Log::debug, *pLog_) << Sm1 << flush;
557 XTP_LOG(Log::debug, *pLog_)
558 << "---------------------------------------" << flush;
559 }
560 XTP_LOG(Log::debug, *pLog_)
561 << TimeStamp() << " Smallest value of dimer overlapmatrix is "
562 << es.eigenvalues()(0) << flush;
563 return J_ortho;
564}
565
566template <class BSE_OPERATOR>
567std::array<Eigen::MatrixXd, 2> BSECoupling::ProjectExcitons(
568 Eigen::MatrixXd& FE_AB, Eigen::MatrixXd& CTStates, BSE_OPERATOR H) const {
569
570 Eigen::MatrixXd projection = OrthogonalizeCTs(FE_AB, CTStates);
571 Eigen::MatrixXd J_ortho = CalcJ_dimer(H, projection);
572
573 std::array<Eigen::MatrixXd, 2> J;
574
575 XTP_LOG(Log::info, *pLog_)
576 << TimeStamp() << " Running Perturbation algorithm" << flush;
577 J[0] = Perturbation(J_ortho);
578 XTP_LOG(Log::info, *pLog_)
579 << TimeStamp() << " Running Projection algorithm" << flush;
580 J[1] = Fulldiag(J_ortho);
581
582 XTP_LOG(Log::debug, *pLog_)
583 << "---------------------------------------" << flush;
584 XTP_LOG(Log::debug, *pLog_) << "Jeff_pert[Hrt]" << flush;
585 XTP_LOG(Log::debug, *pLog_) << J[0] << flush;
586 XTP_LOG(Log::debug, *pLog_) << "Jeff_diag[Hrt]" << flush;
587 XTP_LOG(Log::debug, *pLog_) << J[1] << flush;
588 XTP_LOG(Log::debug, *pLog_)
589 << "---------------------------------------" << flush;
590
591 return J;
592}
593
594Eigen::MatrixXd BSECoupling::Perturbation(
595 const Eigen::MatrixXd& J_dimer) const {
596 Index bse_exc = levA_ + levB_;
597 Index ct = J_dimer.rows() - bse_exc;
598 Eigen::MatrixXd J_result = J_dimer;
599 if (ct > 0) {
600 Eigen::MatrixXd transformation =
601 Eigen::MatrixXd::Identity(J_dimer.rows(), J_dimer.cols());
602 Eigen::MatrixXd Ct = J_dimer.bottomRightCorner(ct, ct);
603 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(Ct);
604 transformation.bottomRightCorner(ct, ct) = es.eigenvectors();
605 Ct.resize(0, 0);
606
607 XTP_LOG(Log::debug, *pLog_) << "FE state hamiltonian" << flush;
608 XTP_LOG(Log::debug, *pLog_)
609 << J_dimer.topLeftCorner(bse_exc, bse_exc) << flush;
610 if (ct > 0) {
611 XTP_LOG(Log::debug, *pLog_) << "eigenvalues of CT states" << flush;
612 XTP_LOG(Log::debug, *pLog_) << es.eigenvalues().transpose() << flush;
613 }
614
615 J_result = transformation.transpose() * J_dimer * transformation;
616 XTP_LOG(Log::debug, *pLog_)
617 << "---------------------------------------" << flush;
618 XTP_LOG(Log::debug, *pLog_) << " J_ortho_[Hrt] CT-state diag" << flush;
619 XTP_LOG(Log::debug, *pLog_) << J_result << flush;
620 XTP_LOG(Log::debug, *pLog_)
621 << "---------------------------------------" << flush;
622 }
623
624 Eigen::MatrixXd Jmatrix = Eigen::MatrixXd::Zero(bse_exc, bse_exc);
625 for (Index stateA = 0; stateA < levA_; stateA++) {
626 double Ea = J_result(stateA, stateA);
627 for (Index stateB = 0; stateB < levB_; stateB++) {
628 Index stateBd = stateB + levA_;
629 XTP_LOG(Log::info, *pLog_)
630 << TimeStamp() << " Calculating coupling between exciton A"
631 << stateA + 1 << " and exciton B" << stateB + 1 << flush;
632 double J = J_result(stateA, stateBd);
633
634 double Eb = J_result(stateBd, stateBd);
635 for (Index k = bse_exc; k < (bse_exc + ct); k++) {
636 double Eab = J_result(k, k);
637 if (std::abs(Eab - Ea) < 0.001) {
638 XTP_LOG(Log::error, *pLog_)
639 << TimeStamp() << "Energydifference between state A "
640 << stateA + 1 << "and CT state " << k + 1 << " is " << Eab - Ea
641 << "[Hrt]" << flush;
642 }
643 if (std::abs(Eab - Eb) < 0.001) {
644 XTP_LOG(Log::error, *pLog_)
645 << TimeStamp() << "Energydifference between state B "
646 << stateB + 1 << "and CT state " << k + 1 << " is " << Eab - Eb
647 << "[Hrt]" << flush;
648 }
649 J += 0.5 * J_result(k, stateA) * J_result(k, stateBd) *
650 (1 / (Ea - Eab) + 1 / (Eb - Eab)); // Have no clue why 0.5
651 }
652 Jmatrix(stateA, stateBd) = J;
653 Jmatrix(stateBd, stateA) = J;
654 }
655 }
656 return Jmatrix;
657}
658
659Eigen::MatrixXd BSECoupling::Fulldiag(const Eigen::MatrixXd& J_dimer) const {
660 Index bse_exc = levA_ + levB_;
661 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(J_dimer);
662 XTP_LOG(Log::debug, *pLog_)
663 << "---------------------------------------" << flush;
664 XTP_LOG(Log::debug, *pLog_) << "Eigenvectors of J" << flush;
665 XTP_LOG(Log::debug, *pLog_) << es.eigenvectors() << flush;
666 XTP_LOG(Log::debug, *pLog_) << "J_eigenvalues[Hrt]" << flush;
667 XTP_LOG(Log::debug, *pLog_) << es.eigenvalues() << flush;
668 XTP_LOG(Log::debug, *pLog_)
669 << "---------------------------------------" << flush;
670 Eigen::MatrixXd Jmat = Eigen::MatrixXd::Zero(bse_exc, bse_exc);
671 // Calculate projection on subspace for every pair of excitons separately
672 for (Index stateA = 0; stateA < levA_; stateA++) {
673 for (Index stateB = 0; stateB < levB_; stateB++) {
674
675 Index stateBd = stateB + levA_;
676 XTP_LOG(Log::info, *pLog_)
677 << TimeStamp() << " Calculating coupling between exciton A"
678 << stateA + 1 << " and exciton B" << stateB + 1 << flush;
679
680 std::array<int, 2> indexes;
681 std::array<int, 2> signs;
682
683 // Find the eigenstate state, which in an L2 is closed to A or B
684 // respectively
685 es.eigenvectors().row(stateA).cwiseAbs().maxCoeff(&indexes[0]);
686 es.eigenvectors().row(stateBd).cwiseAbs().maxCoeff(&indexes[1]);
687 if (indexes[0] == indexes[1]) {
688 Eigen::RowVectorXd stateamplitudes =
689 es.eigenvectors().row(stateBd).cwiseAbs();
690 stateamplitudes[indexes[1]] = 0.0;
691 stateamplitudes.maxCoeff(&indexes[1]);
692 }
693
694 signs[0] = GetSign(es.eigenvectors()(stateA, indexes[0]));
695 signs[1] = GetSign(es.eigenvectors()(stateBd, indexes[1]));
696
697 XTP_LOG(Log::info, *pLog_)
698 << TimeStamp() << " Order is: [Initial state n->nth eigenvalue]"
699 << flush;
700 XTP_LOG(Log::info, *pLog_) << " A" << stateA + 1 << ":" << stateA + 1
701 << "->" << indexes[0] + 1 << " ";
702 XTP_LOG(Log::info, *pLog_) << " B" << stateB + 1 << ":" << stateBd + 1
703 << "->" << indexes[1] + 1 << " " << flush;
704
705 // setting up transformation matrix Tmat and diagonal matrix Emat for the
706 // eigenvalues;
707 Eigen::Matrix2d Emat = Eigen::Matrix2d::Zero();
708 Eigen::Matrix2d Tmat = Eigen::Matrix2d::Zero();
709 // find the eigenvectors which are most similar to the initial states
710 // row
711 for (Index i = 0; i < 2; i++) {
712 Index k = indexes[i];
713 double sign = signs[i];
714 Tmat(0, i) = sign * es.eigenvectors()(stateA, k);
715 Tmat(1, i) = sign * es.eigenvectors()(stateBd, k);
716 Emat(i, i) = es.eigenvalues()(k);
717 }
718 Tmat.colwise().normalize();
719
720 if (Tmat.determinant() < 0) {
721 XTP_LOG(Log::info, *pLog_)
722 << " Reduced state matrix is not in a right handed basis, "
723 "multiplying second eigenvector by -1 "
724 << flush;
725 Tmat.col(1) *= -1;
726 }
727
728 XTP_LOG(Log::debug, *pLog_)
729 << "---------------------------------------" << flush;
730 XTP_LOG(Log::debug, *pLog_) << " T_" << flush;
731 XTP_LOG(Log::debug, *pLog_) << Tmat << flush;
732
733 Eigen::Matrix2d S_small = Tmat * Tmat.transpose();
734
735 XTP_LOG(Log::debug, *pLog_) << "S_small" << flush;
736 XTP_LOG(Log::debug, *pLog_) << S_small << flush;
737 // orthogonalize that matrix
738
739 Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> ss(S_small);
740 Eigen::Matrix2d sm1 = ss.operatorInverseSqrt();
741 Emat = sm1 * Emat * sm1;
742
743 XTP_LOG(Log::info, *pLog_)
744 << TimeStamp() << " Smallest value of dimer overlapmatrix is "
745 << ss.eigenvalues()(0) << flush;
746
747 XTP_LOG(Log::debug, *pLog_) << "S-1/2" << flush;
748 XTP_LOG(Log::debug, *pLog_) << sm1 << flush;
749 XTP_LOG(Log::debug, *pLog_) << "E_ortho" << flush;
750 XTP_LOG(Log::debug, *pLog_) << Emat << flush;
751
752 Tmat = Tmat * sm1;
753
754 XTP_LOG(Log::debug, *pLog_) << "T_ortho" << flush;
755 XTP_LOG(Log::debug, *pLog_) << Tmat << flush;
756 XTP_LOG(Log::debug, *pLog_)
757 << "---------------------------------------" << flush;
758
759 Eigen::Matrix2d J_small = Tmat * Emat * Tmat.transpose();
760 XTP_LOG(Log::debug, *pLog_) << "T_ortho*E_ortho*T_ortho^T" << flush;
761 XTP_LOG(Log::debug, *pLog_) << J_small << flush;
762
763 Jmat(stateA, stateBd) = J_small(0, 1);
764 Jmat(stateBd, stateA) = J_small(1, 0);
765 }
766 }
767
768 return Jmat;
769}
770
771} // namespace xtp
772} // namespace votca
bse_operator.h
votca::tools::EigenSystem::eigenvalues
const Eigen::VectorXd & eigenvalues() const
Definition eigensystem.h:30
votca::tools::EigenSystem::eigenvectors
const Eigen::MatrixXd & eigenvectors() const
Definition eigensystem.h:33
votca::tools::Property
class to manage program options with xml serialization functionality
Definition property.h:55
votca::tools::Property::add
Property & add(const std::string &key, const std::string &value)
add a new property to structure
Definition property.cc:108
votca::tools::Property::get
Property & get(const std::string &key)
get existing property
Definition property.cc:79
votca::tools::Property::as
T as() const
return value as type
Definition property.h:283
votca::tools::Property::setAttribute
void setAttribute(const std::string &attribute, const T &value)
set an attribute
Definition property.h:314
votca::xtp::AOBasis
Container to hold Basisfunctions for all atoms.
Definition aobasis.h:42
votca::xtp::AOBasis::AOBasisSize
Index AOBasisSize() const
Definition aobasis.h:46
votca::xtp::BSECoupling::getSingletCouplingElement
double getSingletCouplingElement(Index levelA, Index levelB, Index methodindex) const
Definition bsecoupling.cc:122
votca::xtp::BSECoupling::getTripletCouplingElement
double getTripletCouplingElement(Index levelA, Index levelB, Index methodindex) const
Definition bsecoupling.cc:128
votca::xtp::BSECoupling::Identify
std::string Identify() const
Definition bsecoupling.h:43
votca::xtp::BSECoupling::CalculateCouplings
void CalculateCouplings(const Orbitals &orbitalsA, const Orbitals &orbitalsB, const Orbitals &orbitalsAB) override
evaluates electronic couplings
Definition bsecoupling.cc:220
votca::xtp::BSECoupling::Addoutput
void Addoutput(tools::Property &type_summary, const Orbitals &orbitalsA, const Orbitals &orbitalsB) const override
Definition bsecoupling.cc:88
votca::xtp::BSECoupling::WriteToProperty
void WriteToProperty(tools::Property &summary, const QMState &stateA, const QMState &stateB) const
Definition bsecoupling.cc:66
votca::xtp::BSECoupling::JAB_triplet
std::array< Eigen::MatrixXd, 2 > JAB_triplet
Definition bsecoupling.h:95
votca::xtp::BSECoupling::Perturbation
Eigen::MatrixXd Perturbation(const Eigen::MatrixXd &J_dimer) const
Definition bsecoupling.cc:594
votca::xtp::BSECoupling::Initialize
void Initialize(tools::Property &options) override
Definition bsecoupling.cc:38
votca::xtp::BSECoupling::JAB_singlet
std::array< Eigen::MatrixXd, 2 > JAB_singlet
Definition bsecoupling.h:94
votca::xtp::BSECoupling::ProjectFrenkelExcitons
Eigen::MatrixXd ProjectFrenkelExcitons(const Eigen::MatrixXd &BSE_Coeffs, const Eigen::MatrixXd &X_AB, Index bseX_vtotal, Index bseX_ctotal, Index bseAB_vtotal, Index bseAB_ctotal) const
Definition bsecoupling.cc:190
votca::xtp::BSECoupling::SetupCTStates
Eigen::MatrixXd SetupCTStates(Index bseA_vtotal, Index bseB_vtotal, Index bseAB_vtotal, Index bseAB_ctotal, const Eigen::MatrixXd &A_AB, const Eigen::MatrixXd &B_AB) const
Definition bsecoupling.cc:134
votca::xtp::BSECoupling::ProjectExcitons
std::array< Eigen::MatrixXd, 2 > ProjectExcitons(Eigen::MatrixXd &FE_AB, Eigen::MatrixXd &CTStates, BSE_OPERATOR H) const
Definition bsecoupling.cc:567
votca::xtp::BSECoupling::Fulldiag
Eigen::MatrixXd Fulldiag(const Eigen::MatrixXd &J_dimer) const
Definition bsecoupling.cc:659
votca::xtp::BSECoupling::CalcJ_dimer
Eigen::MatrixXd CalcJ_dimer(BSE_OPERATOR &H, Eigen::MatrixXd &projection) const
Definition bsecoupling.cc:511
votca::xtp::BSECoupling::OrthogonalizeCTs
Eigen::MatrixXd OrthogonalizeCTs(Eigen::MatrixXd &FE_AB, Eigen::MatrixXd &CTStates) const
Definition bsecoupling.cc:471
votca::xtp::BSE_OPERATOR
Definition bse_operator.h:41
votca::xtp::BSE::getTripletOperator_TDA
TripletOperator_TDA getTripletOperator_TDA() const
Definition bse.cc:173
votca::xtp::BSE::getSingletOperator_TDA
SingletOperator_TDA getSingletOperator_TDA() const
Definition bse.cc:166
votca::xtp::BSE::configure
void configure(const options &opt, const Eigen::VectorXd &RPAEnergies, const Eigen::MatrixXd &Hqp_in)
Definition bse.cc:43
votca::xtp::CouplingBase::CalculateOverlapMatrix
Eigen::MatrixXd CalculateOverlapMatrix(const Orbitals &orbitalsAB) const
Definition couplingbase.cc:27
votca::xtp::CouplingBase::CheckAtomCoordinates
void CheckAtomCoordinates(const Orbitals &orbitalsA, const Orbitals &orbitalsB, const Orbitals &orbitalsAB) const
Definition couplingbase.cc:35
votca::xtp::Orbitals
container for molecular orbitals
Definition orbitals.h:46
votca::xtp::Orbitals::getBSEvmax
Index getBSEvmax() const
Definition orbitals.h:286
votca::xtp::Orbitals::BSETriplets
const tools::EigenSystem & BSETriplets() const
Definition orbitals.h:324
votca::xtp::Orbitals::GetFlagUseHqpOffdiag
bool GetFlagUseHqpOffdiag() const
Definition orbitals.h:408
votca::xtp::Orbitals::getBSEcmin
Index getBSEcmin() const
Definition orbitals.h:288
votca::xtp::Orbitals::QPdiag
const tools::EigenSystem & QPdiag() const
Definition orbitals.h:317
votca::xtp::Orbitals::getHomo
Index getHomo() const
Definition orbitals.h:68
votca::xtp::Orbitals::BSESinglets
const tools::EigenSystem & BSESinglets() const
Definition orbitals.h:334
votca::xtp::Orbitals::getBSEcmax
Index getBSEcmax() const
Definition orbitals.h:290
votca::xtp::Orbitals::getAuxBasis
const AOBasis & getAuxBasis() const
Definition orbitals.h:218
votca::xtp::Orbitals::getBSEvmin
Index getBSEvmin() const
Definition orbitals.h:284
votca::xtp::Orbitals::getRPAmax
Index getRPAmax() const
Definition orbitals.h:264
votca::xtp::Orbitals::getBasisSetSize
Index getBasisSetSize() const
Definition orbitals.h:64
votca::xtp::Orbitals::RPAInputEnergies
const Eigen::VectorXd & RPAInputEnergies() const
Definition orbitals.h:299
votca::xtp::Orbitals::getGWAmin
Index getGWAmin() const
Definition orbitals.h:249
votca::xtp::Orbitals::MOs
const tools::EigenSystem & MOs() const
Definition orbitals.h:122
votca::xtp::Orbitals::getGWAmax
Index getGWAmax() const
Definition orbitals.h:251
votca::xtp::Orbitals::getDFTbasisName
const std::string & getDFTbasisName() const
Definition orbitals.h:202
votca::xtp::Orbitals::getRPAmin
Index getRPAmin() const
Definition orbitals.h:262
votca::xtp::Orbitals::getDftBasis
const AOBasis & getDftBasis() const
Definition orbitals.h:208
votca::xtp::QMStateType::ToLongString
std::string ToLongString() const
Definition qmstate.cc:69
votca::xtp::QMState
Identifier for QMstates. Strings like S1 are converted into enum +zero indexed int.
Definition qmstate.h:132
votca::xtp::QMState::ToString
std::string ToString() const
Definition qmstate.cc:146
votca::xtp::QMState::StateIdx
Index StateIdx() const
Definition qmstate.h:154
votca::xtp::QMState::Type
const QMStateType & Type() const
Definition qmstate.h:151
votca::xtp::TCMatrix_gwbse
Definition threecenter.h:69
votca::xtp::TCMatrix_gwbse::Fill
void Fill(const AOBasis &auxbasis, const AOBasis &dftbasis, const Eigen::MatrixXd &dft_orbitals)
Definition threecenter.cc:72
votca::xtp::TCMatrix_gwbse::Initialize
void Initialize(Index basissize, Index mmin, Index mmax, Index nmin, Index nmax)
Definition threecenter.cc:29
votca::xtp::TimeStamp
Timestamp returns the current time as a string Example: cout << TimeStamp()
Definition logger.h:224
XTP_LOG
#define XTP_LOG(level, log)
Definition logger.h:40
std
STL namespace.
votca::tools::conv::hrt2ev
const double hrt2ev
Definition constants.h:53
votca::xtp::OPENMP::getMaxThreads
Index getMaxThreads()
Definition eigen.h:128
votca::xtp::GetSign
int GetSign(double value)
Definition bsecoupling.cc:211
votca
base class for all analysis tools
Definition basebead.h:33
votca::Index
Eigen::Index Index
Definition types.h:26
Generated by doxygen 1.9.7