votca 2024.2-dev
Loading...
Searching...
No Matches
ERIs.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// Local VOTCA includes
21#include "votca/xtp/ERIs.h"
22#include "votca/xtp/aobasis.h"
24namespace votca {
25namespace xtp {
26
27void ERIs::Initialize(const AOBasis& dftbasis, const AOBasis& auxbasis) {
28 threecenter_.Fill(auxbasis, dftbasis);
29 return;
30}
31
32void ERIs::Initialize_4c(const AOBasis& dftbasis) {
33
34 basis_ = dftbasis.GenerateLibintBasis();
35 shellpairs_ = dftbasis.ComputeShellPairs();
37 maxnprim_ = dftbasis.getMaxNprim();
38 maxL_ = dftbasis.getMaxL();
39
41
43 return;
44}
45
46std::vector<std::vector<libint2::ShellPair>> ERIs::ComputeShellPairData(
47 const std::vector<libint2::Shell>& basis,
48 const std::vector<std::vector<Index>>& shellpairs) const {
49 std::vector<std::vector<libint2::ShellPair>> shellpairdata(basis.size());
50 const double ln_max_engine_precision =
51 std::log(std::numeric_limits<double>::epsilon() * 1e-10);
52
53#pragma omp parallel for schedule(dynamic)
54 for (Index s1 = 0; s1 < Index(shellpairs.size()); s1++) {
55 for (Index s2 : shellpairs[s1]) {
56 shellpairdata[s1].emplace_back(
57 libint2::ShellPair(basis[s1], basis[s2], ln_max_engine_precision));
58 }
59 }
60 return shellpairdata;
61}
62
63Eigen::MatrixXd ERIs::ComputeShellBlockNorm(const Eigen::MatrixXd& dmat) const {
64 Eigen::MatrixXd result =
65 Eigen::MatrixXd::Zero(starts_.size(), starts_.size());
66#pragma omp parallel for schedule(dynamic)
67 for (Index s1 = 0l; s1 < Index(basis_.size()); ++s1) {
68 Index bf1 = starts_[s1];
69 Index n1 = basis_[s1].size();
70 for (Index s2 = 0l; s2 <= s1; ++s2) {
71 Index bf2 = starts_[s2];
72 Index n2 = basis_[s2].size();
73
74 result(s2, s1) = dmat.block(bf2, bf1, n2, n1).cwiseAbs().maxCoeff();
75 }
76 }
77 return result.selfadjointView<Eigen::Upper>();
78}
79
80Eigen::MatrixXd ERIs::CalculateERIs_3c(const Eigen::MatrixXd& DMAT) const {
81 assert(threecenter_.size() > 0 &&
82 "Please call Initialize before running this");
83 Eigen::MatrixXd ERIs2 = Eigen::MatrixXd::Zero(DMAT.rows(), DMAT.cols());
84 Symmetric_Matrix dmat_sym = Symmetric_Matrix(DMAT);
85#pragma omp parallel for schedule(guided) reduction(+ : ERIs2)
86 for (Index i = 0; i < threecenter_.size(); i++) {
87 const Symmetric_Matrix& threecenter = threecenter_[i];
88 // Trace over prod::DMAT,I(l)=componentwise product over
89 const double factor = threecenter.TraceofProd(dmat_sym);
90 Eigen::SelfAdjointView<Eigen::MatrixXd, Eigen::Upper> m =
91 ERIs2.selfadjointView<Eigen::Upper>();
92 threecenter.AddtoEigenUpperMatrix(m, factor);
93 }
94
95 return ERIs2.selfadjointView<Eigen::Upper>();
96}
97
98Eigen::MatrixXd ERIs::CalculateEXX_dmat(const Eigen::MatrixXd& DMAT) const {
99 assert(threecenter_.size() > 0 &&
100 "Please call Initialize before running this");
101 Eigen::MatrixXd EXX = Eigen::MatrixXd::Zero(DMAT.rows(), DMAT.cols());
102
103#pragma omp parallel for schedule(guided) reduction(+ : EXX)
104 for (Index i = 0; i < threecenter_.size(); i++) {
105 const Eigen::MatrixXd threecenter = threecenter_[i].UpperMatrix();
106 EXX -= threecenter.selfadjointView<Eigen::Upper>() * DMAT *
107 threecenter.selfadjointView<Eigen::Upper>();
108 }
109 return EXX;
110}
111
112Eigen::MatrixXd ERIs::CalculateEXX_mos(const Eigen::MatrixXd& occMos) const {
113 assert(threecenter_.size() > 0 &&
114 "Please call Initialize before running this");
115 Eigen::MatrixXd EXX = Eigen::MatrixXd::Zero(occMos.rows(), occMos.rows());
116
117#pragma omp parallel for schedule(guided) reduction(+ : EXX)
118 for (Index i = 0; i < threecenter_.size(); i++) {
119 const Eigen::MatrixXd TCxMOs_T =
120 occMos.transpose() *
121 threecenter_[i].UpperMatrix().selfadjointView<Eigen::Upper>();
122 EXX -= TCxMOs_T.transpose() * TCxMOs_T;
123 }
124 return 2 * EXX;
125}
126
127} // namespace xtp
128} // namespace votca
Container to hold Basisfunctions for all atoms.
Definition aobasis.h:42
Index getMaxL() const
Definition aobasis.cc:29
std::vector< Index > getMapToBasisFunctions() const
Definition aobasis.cc:45
Index getMaxNprim() const
Definition aobasis.cc:37
std::vector< std::vector< Index > > ComputeShellPairs(double threshold=1e-20) const
std::vector< libint2::Shell > GenerateLibintBasis() const
Definition aobasis.cc:119
Eigen::MatrixXd CalculateERIs_3c(const Eigen::MatrixXd &DMAT) const
Definition ERIs.cc:80
std::vector< std::vector< libint2::ShellPair > > shellpairdata_
Definition ERIs.h:78
std::vector< libint2::Shell > basis_
Definition ERIs.h:74
Eigen::MatrixXd ComputeSchwarzShells(const AOBasis &dftbasis) const
std::vector< std::vector< Index > > shellpairs_
Definition ERIs.h:77
void Initialize_4c(const AOBasis &dftbasis)
Definition ERIs.cc:32
Index maxL_
Definition ERIs.h:80
std::vector< std::vector< libint2::ShellPair > > ComputeShellPairData(const std::vector< libint2::Shell > &basis, const std::vector< std::vector< Index > > &shellpairs) const
Definition ERIs.cc:46
void Initialize(const AOBasis &dftbasis, const AOBasis &auxbasis)
Definition ERIs.cc:27
Eigen::MatrixXd schwarzscreen_
Definition ERIs.h:98
Eigen::MatrixXd ComputeShellBlockNorm(const Eigen::MatrixXd &dmat) const
Definition ERIs.cc:63
std::vector< Index > starts_
Definition ERIs.h:75
Eigen::MatrixXd CalculateEXX_dmat(const Eigen::MatrixXd &DMAT) const
Definition ERIs.cc:98
TCMatrix_dft threecenter_
Definition ERIs.h:96
Index maxnprim_
Definition ERIs.h:79
Eigen::MatrixXd CalculateEXX_mos(const Eigen::MatrixXd &occMos) const
Definition ERIs.cc:112
void AddtoEigenUpperMatrix(Eigen::SelfAdjointView< Eigen::MatrixXd, Eigen::Upper > &upper, double factor=1.0) const
double TraceofProd(const Symmetric_Matrix &a) const
void Fill(const AOBasis &auxbasis, const AOBasis &dftbasis)
base class for all analysis tools
Definition basebead.h:33
Eigen::Index Index
Definition types.h:26