Program Listing for File ri_utils.hpp

↰ Return to documentation for file (src/three_electron/ri_utils.hpp)

//
// Created by Zack Williams on 25/03/2024.
//

#ifndef RI_UTILS_HPP
#define RI_UTILS_HPP

#include "utils/linalg.hpp"

namespace uw12::three_el::ri {

struct ABSProjectors {
  linalg::Mat s_inv_ri_ri;
  linalg::Mat s_inv_ri_ao;
  linalg::Mat s_inv_ao_ri;
  linalg::Mat s_inv_ao_ao;
};

inline ABSProjectors calculate_abs_projectors(
    const linalg::Mat& overlap_matrix,
    const size_t n_ao,
    const size_t n_ri,
    const double eigen_ld_threshold = 1e-8,
    const double linear_dependency_threshold = 1e-6
) {
  using linalg::head_rows;
  using linalg::tail_rows;
  using linalg::transpose;

  const auto [s_vals, s_vecs] = linalg::eigen_decomposition(
      overlap_matrix, linear_dependency_threshold, eigen_ld_threshold
  );

  if (linalg::n_rows(s_vecs) != n_ao + n_ri) {
    throw std::runtime_error("Total number of basis functions is inconsistent");
  }

  assert(linalg::n_cols(s_vecs) == linalg::n_elem(s_vals));

  const linalg::Mat inv = linalg::inv_sym_pd(linalg::diagmat(s_vals));

  return {
      tail_rows(s_vecs, n_ri) * inv * transpose(tail_rows(s_vecs, n_ri)),
      tail_rows(s_vecs, n_ri) * inv * transpose(head_rows(s_vecs, n_ao)),
      head_rows(s_vecs, n_ao) * inv * transpose(tail_rows(s_vecs, n_ri)),
      head_rows(s_vecs, n_ao) * inv * transpose(head_rows(s_vecs, n_ao))
  };
}

}  // namespace uw12::three_el::ri
#endif  // RI_UTILS_HPP