/*
* Copyright 2009-2020 The VOTCA Development Team
* (http://www.votca.org)
*
* Licensed under the Apache License, Version 2.0 (the "License")
*
* You may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*/
#include <votca/tools/constants.h>
#include "votca/xtp/ImaginaryAxisIntegration.h"
#include "votca/xtp/quadrature_factory.h"
#include "votca/xtp/threecenter.h"
namespace votca {
namespace xtp {
ImaginaryAxisIntegration::ImaginaryAxisIntegration(
const Eigen::VectorXd& energies, const TCMatrix_gwbse& Mmn)
: _energies(energies), _Mmn(Mmn) {}
void ImaginaryAxisIntegration::configure(
options opt, const RPA& rpa, const Eigen::MatrixXd& kDielMxInv_zero) {
_opt = opt;
QuadratureFactory::RegisterAll();
_gq = std::unique_ptr<GaussianQuadratureBase>(
Quadratures().Create(_opt.quadrature_scheme));
_gq->configure(_opt.order);
CalcDielInvVector(rpa, kDielMxInv_zero);
}
// This function calculates and stores inverses of the microscopic dielectric
// matrix in a matrix vector
void ImaginaryAxisIntegration::CalcDielInvVector(
const RPA& rpa, const Eigen::MatrixXd& kDielMxInv_zero) {
_dielinv_matrices_r.resize(_gq->Order());
for (Index j = 0; j < _gq->Order(); j++) {
double newpoint = _gq->ScaledPoint(j);
Eigen::MatrixXd eps_inv_j = rpa.calculate_epsilon_i(newpoint).inverse();
eps_inv_j.diagonal().array() -= 1.0;
_dielinv_matrices_r[j] =
-eps_inv_j +
kDielMxInv_zero * std::exp(-std::pow(_opt.alpha * newpoint, 2));
}
}
class FunctionEvaluation {
public:
FunctionEvaluation(const Eigen::MatrixXd& Imx, const Eigen::ArrayXcd& DeltaE,
const std::vector<Eigen::MatrixXd>& dielinv_matrices_r)
: _Imx(Imx), _DeltaE(DeltaE), _dielinv_matrices_r(dielinv_matrices_r){};
double operator()(Index j, double point, bool symmetry) const {
Eigen::VectorXcd denominator;
const std::complex<double> cpoint(0.0, point);
if (symmetry) {
denominator =
(_DeltaE + cpoint).cwiseInverse() + (_DeltaE - cpoint).cwiseInverse();
} else {
denominator = (_DeltaE + cpoint).cwiseInverse();
}
return 0.5 / tools::conv::Pi *
((_Imx * (_dielinv_matrices_r[j].conjugate()))
.cwiseProduct(denominator.asDiagonal() * _Imx))
.sum()
.real();
}
private:
const Eigen::MatrixXd& _Imx;
const Eigen::ArrayXcd& _DeltaE;
const std::vector<Eigen::MatrixXd>& _dielinv_matrices_r;
};
double ImaginaryAxisIntegration::SigmaGQDiag(double frequency, Index gw_level,
double eta) const {
Index lumo = _opt.homo + 1;
const Index occ = lumo - _opt.rpamin;
const Index unocc = _opt.rpamax - _opt.homo;
Index gw_level_offset = gw_level + _opt.qpmin - _opt.rpamin;
const Eigen::MatrixXd& Imx = _Mmn[gw_level_offset];
Eigen::ArrayXcd DeltaE = frequency - _energies.array();
DeltaE.imag().head(occ) = eta;
DeltaE.imag().tail(unocc) = -eta;
FunctionEvaluation f = FunctionEvaluation(Imx, DeltaE, _dielinv_matrices_r);
return _gq->Integrate(f);
}
} // namespace xtp
} // namespace votca