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 #include "sigma_cda.h"
21 #include "votca/xtp/gw.h"
22 #include <votca/tools/constants.h>
23
24 namespace votca {
25 namespace xtp {
26
27 // Prepares the zero and imaginary frequency kappa matrices with
28 // kappa(omega) = epsilon^-1(omega) - 1 needed in numerical
29 // integration and for the Gaussian tail
PrepareScreening()30 void Sigma_CDA::PrepareScreening() {
31 ImaginaryAxisIntegration::options opt;
32 opt.homo = opt_.homo;
33 opt.order = opt_.order;
34 opt.qptotal = qptotal_;
35 opt.qpmin = opt_.qpmin;
36 opt.rpamax = opt_.rpamax;
37 opt.rpamin = opt_.rpamin;
38 opt.alpha = opt_.alpha;
39 opt.quadrature_scheme = opt_.quadrature_scheme;
40 // prepare the zero frequency inverse for Gaussian tail
41 kDielMxInv_zero_ =
42 rpa_.calculate_epsilon_r(std::complex<double>(0.0, 0.0)).inverse();
43 kDielMxInv_zero_.diagonal().array() -= 1.0;
44 gq_.configure(opt, rpa_, kDielMxInv_zero_);
45 }
46
47 // This function is used in the calculation of the residues and
48 // calculates the real part of the dielectric function for a complex
49 // frequency of the kind omega = delta + i*eta. Instead of explicit
50 // inversion and multiplication with and Imx vector, a linear system
51 // is solved.
CalcDiagContribution(const Eigen::MatrixXd::ConstRowXpr & Imx_row,double delta,double eta) const52 double Sigma_CDA::CalcDiagContribution(
53 const Eigen::MatrixXd::ConstRowXpr& Imx_row, double delta,
54 double eta) const {
55 std::complex<double> delta_eta(delta, eta);
56
57 Eigen::MatrixXd DielMxInv = rpa_.calculate_epsilon_r(delta_eta);
58 Eigen::VectorXd x =
59 DielMxInv.partialPivLu().solve(Imx_row.transpose()) - Imx_row.transpose();
60 return x.dot(Imx_row.transpose());
61 }
62
63 // Step-function prefactor for the residues
CalcResiduePrefactor(double e_f,double e_m,double frequency) const64 double Sigma_CDA::CalcResiduePrefactor(double e_f, double e_m,
65 double frequency) const {
66 double factor = 0.0;
67 double tolerance = 1e-10;
68 if (e_f < e_m && e_m < frequency) {
69 factor = 1.0;
70 } else if (e_f > e_m && e_m > frequency) {
71 factor = -1.0;
72 } else if (std::abs(e_m - frequency) < tolerance && e_f > e_m) {
73 factor = -0.5;
74 } else if (std::abs(e_m - frequency) < tolerance && e_f < e_m) {
75 factor = 0.5;
76 }
77 return factor;
78 }
79
80 // Calculates the contribution of residues to the correlation
81 // part of the self-energy of a fixed gw_level
CalcResidueContribution(double frequency,Index gw_level) const82 double Sigma_CDA::CalcResidueContribution(double frequency,
83 Index gw_level) const {
84
85 const Eigen::VectorXd& rpa_energies = rpa_.getRPAInputEnergies();
86 Index rpatotal = rpa_energies.size();
87 Index gw_level_offset = gw_level + opt_.qpmin - opt_.rpamin;
88
89 double sigma_c = 0.0;
90 double sigma_c_tail = 0.0;
91 Index homo = opt_.homo - opt_.rpamin;
92 Index lumo = homo + 1;
93 double fermi_rpa = (rpa_energies(lumo) + rpa_energies(homo)) / 2.0;
94 const Eigen::MatrixXd& Imx = Mmn_[gw_level_offset];
95
96 for (Index i = 0; i < rpatotal; ++i) {
97 double delta = rpa_energies(i) - frequency;
98 double abs_delta = std::abs(delta);
99 double factor = CalcResiduePrefactor(fermi_rpa, rpa_energies(i), frequency);
100
101 // Only considering the terms with a abs(prefactor) > 0.
102 // The prefactor can be 1,-1,0.5,-0.5 or 0. We avoid calculating the
103 // diagonal contribution if the prefactor is 0. We want to calculate it for
104 // all the other cases.
105 if (std::abs(factor) > 1e-10) {
106 sigma_c +=
107 factor * CalcDiagContribution(Imx.row(i), abs_delta, rpa_.getEta());
108 }
109 // adds the contribution from the Gaussian tail
110 if (abs_delta > 1e-10) {
111 sigma_c_tail +=
112 CalcDiagContributionValue_tail(Imx.row(i), delta, opt_.alpha);
113 }
114 }
115 return sigma_c + sigma_c_tail;
116 }
117
118 // Calculates the correlation part of the self-energy for a fixed
119 // gw_level and frequence by evaluating the numerical integration
120 // and residue contributions
CalcCorrelationDiagElement(Index gw_level,double frequency) const121 double Sigma_CDA::CalcCorrelationDiagElement(Index gw_level,
122 double frequency) const {
123
124 double sigma_c_residue = CalcResidueContribution(frequency, gw_level);
125 double sigma_c_integral = gq_.SigmaGQDiag(frequency, gw_level, rpa_.getEta());
126 return sigma_c_residue + sigma_c_integral;
127 }
128
129 // Calculates the contribuion of the tail correction to the
130 // residue term
CalcDiagContributionValue_tail(const Eigen::MatrixXd::ConstRowXpr & Imx_row,double delta,double alpha) const131 double Sigma_CDA::CalcDiagContributionValue_tail(
132 const Eigen::MatrixXd::ConstRowXpr& Imx_row, double delta,
133 double alpha) const {
134
135 double erfc_factor = 0.5 * std::copysign(1.0, delta) *
136 std::exp(std::pow(alpha * delta, 2)) *
137 std::erfc(std::abs(alpha * delta));
138
139 double value = (Imx_row * kDielMxInv_zero_).dot(Imx_row);
140 return value * erfc_factor;
141 }
142
143 } // namespace xtp
144 } // namespace votca
145