24#include <libecpint/ecpint.hpp>
25#include <libecpint/gshell.hpp>
26#include <libecpint/mathutil.hpp>
27#include <libint2/solidharmonics.h>
33 const libecpint::GaussianShell& shell) {
35 <<
" L:" <<
Index(shell.am()) <<
" Func:" << shell.nprimitive() <<
"\n";
36 for (
int i = 0; i < shell.nprimitive(); i++) {
37 out <<
" Gaussian Decay: " << shell.exp(i);
38 out <<
" Contractions: " << shell.coef(i) <<
"\n";
44 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
50 std::vector<libecpint::GaussianShell> basis;
51 std::vector<Index> cartesian_size;
52 std::vector<Index> spherical_size;
53 for (
const auto& shell : aobasis) {
54 libecpint::GaussianShell s(
55 {shell.getPos().x(), shell.getPos().y(), shell.getPos().z()},
57 spherical_size.push_back(shell.getNumFunc());
58 cartesian_size.push_back(shell.getCartesianNumFunc());
61 libint2::Shell s_libint = shell.LibintShell();
62 for (
Index i = 0; i <
Index(s_libint.nprim()); i++) {
63 s.addPrim(s_libint.alpha[i], s_libint.contr[0].coeff[i]);
69 std::vector<libecpint::ECPIntegral> engines(
71 libecpint::ECPIntegral(
int(aobasis.
getMaxL()),
int(ecp.
getMaxL()), 0));
72#pragma omp parallel for schedule(guided)
75 libecpint::ECPIntegral& engine = engines[thread_id];
76 Index bf1 = shell2bf[s1];
77 Index n1 = spherical_size[s1];
78 Index c1 = cartesian_size[s1];
79 for (
Index s2 = 0; s2 <= s1; ++s2) {
80 Index bf2 = shell2bf[s2];
81 Index n2 = spherical_size[s2];
82 Index c2 = cartesian_size[s2];
84 MatrixLibInt cartesian_result = MatrixLibInt::Zero(c1, c2);
85 for (
auto& ecppotential : ecp) {
86 libecpint::TwoIndex<double> results;
87 engine.compute_shell_pair(ecppotential, basis[s1], basis[s2], results);
89 Eigen::Map<MatrixLibInt>(results.data.data(), c1, c2);
92 MatrixLibInt spherical_result = MatrixLibInt::Zero(n1, n2);
93 libint2::solidharmonics::tform<double>(basis[s1].l, basis[s2].l,
94 cartesian_result.data(),
95 spherical_result.data());
99 aopotential_.block(bf2, bf1, n2, n1) = spherical_result.transpose();
Container to hold Basisfunctions for all atoms.
Index AOBasisSize() const
std::vector< Index > getMapToBasisFunctions() const
Index getNumofShells() const
void FillPotential(const AOBasis &aobasis, const ECPAOBasis &ecp)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > aopotential_
Container to hold ECPs for all atoms.
std::string EnumToString(L l)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixLibInt
base class for all analysis tools